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I.  Abstract 

The'stucty  presents  a  method,  oased  on  the  flux  vector  splitting  approach, 
to  the  problem  of  unsteady  two  dimensional  inviscid  transonic  flows,  with 
emphasis  on  the  numerical  determination  of  the  shock  position,  through 
nozzles  with  varying  back  pressure.  The  methodology/ governed  by  'the 
Euler  equations’  is  first  explained  for  one  and  two  dimensional  steady 
state  applications,  and  the  accuracy  of  the  results  is  validated  by 
comparison  with  exact  (one  dimension)  and  numerical  (two  dimensions) 
solutions. 

Tne  model  is  then  applied  to  the  problem  of  fluctuating  back  pressure  in 
quasi  one  dimensional  and  two  dimensional  flows.  The  one  dimensional 
results  are  validated  by  comparison  with  a  small  perturbation  analytical 
unsteady  solution,  whereafter  sample  cases  are  performed  with  the 
objective  to  understand  fundamental  aspects  of  unsteady  flows, 
it  is  concluded  that  both  the  amplitude  and  frequency  of  the  imposed 
fluctuating  exit  pressure  are  important  parameters  for  ^location  of  the 
unsteady  shock,  it  is  aiso  shown  thar  tne  average  unsteady  shock  position 
is  not  identical  with  the  steady  state  position,  and  that  the  unsteady 
shock  may,  under  certain  circumstances,  propagate  upstream  into  the 
subsonic  flow  domain.  The  pressure  jump  over  the  shock,  as  well  as  the 
unsteady  post-shock  pressure,  is  different  for  identical  shock-positions 
during  tne  cycle  of  fluctuation,  which  implies  that  an  unsteady  shock 
movement,  imposed  by  oscillating  back  pressure,  may  introduce  a 
significant  unsteady  lift  and  moment.  This  may  be  of  importance  for 
flutter  predictions. 

it  is  also  noted  that,  although  the  some  velocity  is  obtained  in  tne  throat 
of  steady  state  quasi  one  dimensional  flow,  this  is  not  necessarily  true 
for  the  unsteady  solution.  During  part  of  the  period  with  fluctuating  back 
pressure,  the  flew  velocity  may  be  subsonic  at  the  threat,  and  still  reach 
a  supersonic  value  later  in  the  nozzle.  This  phenomena  depends  on  tne. 
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II.  Brief  Information  to  the  Reader. 


The  present  feport  consists  of  four  main  parts,  each  of  which  can  quickly 
be  scanned  in  dependence  of  the  readers  interest.  The  following  gives  a 
brief  overview  of  the  different  parts. 

A  general  introduction  to  blade  flutter  in  axial  flow  turbomachines  is 
given  in  section  i.i,  followed  by  an  introduction  to  existing  numerical 
methods  for  attacking  this  problem.  It  is  not  necessary  to  go  into  this 
section  in  order  to  judge  the  results  presented  in  the  later  sections,  it 
may  nowever  be  of  interest  if  the  reader  would  like  to  judge  the 
implications  of  tne  results  in  the  broade  context  of  blade  flutter. 

The  numerically  method  employed  in  the  present  study  is  briefly  presented 
in  sections  3.1  to  3 .4  and  4.1  to  4.5.  Again,  the  results  can  be  judged 
essentially  without  reading  these  sections. 

The  validation  of  the  method  and  a  few  sample  cases  are  presented  in 
sections  3.5  (quasi  one  dimensional  steady  state),  4.6  (two  dimensional 
steady  state),  5.1  to  5.2  (quasi  one  dimensional  unsteady)  and  5.3  (two 
dimensional  unsteady).  These  sections  are  obviously  the  most  important 
ones. 

A  few  details  about  the  method  are  given  in  the  appendices,  together  with 
indications  of  how  to  run  the  program  on  the  IBM  main  frame  at  the  Naval 
Postgraduate  School,  Monterey. 

it  should  also  be  stated  that  the  report  treats  "Work  i  progress"  and  is 
largely  meant  to  be  a  reference  also  for  the  author.  Thus,  several  details 
that  wouid  normally  not  be  found  in  a  final  report  may  be  present  here.  On 
the  other  nano,  certain  aspects  of  a  final  report  in  tne  usual  sense  may 
also  be  missing. 
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LL  General  introduction  to  the  Blade  Flutter  Problem. 

One  of  the  many  serious  problems  encountered  during  the  development  of 
modern  Jet  engines  and  Industrial  turbomachines  is  the  prediction  of  flow 
induced  vibrations  which  may  lead  to  an  eventual  machine  failure.  Due  to 
the  interaction  of  several  complicated  phenomena  in  turbomachines,  this 
problem  is  still  today  far  from  being  resolved  even  though  a  significant 
amount  of  worx  has  been  performed  in  the  domain  of  aeroe las ti city,  wmcn 
concerns  the  vibration  of  a  deformable  structure  in  a  fluid,  during  the  two 
last  decenma. 

Among  all  the  complicated  phenomena,  the  ones  concerning  the  aeroelastic 
stability  of  the  turbomachine  blades  are  part  of  the  most  vital.  Such 
stability  investigations  will  be  even  more  significant  in  the  future  as 
blade  vibration  predictions  become  increasingly  important  as  part  of  a 
trend  towards  higher  unit  powers,  often  without  increase  in  the  number  of 
stages  in  the  machine.  This  trend  implies,  as  examples: 
longer  and  more  slender  blades 
higher  blade  loadings 
higher  flow  velocities 
higher  inlet  flow  temperatures  in  turbines 
higher  mass  flow  rates. 

The  excitation  of  blade  vibrations  may  be  either  forced  or  selfsustained. 
in  the  first  group  (forced  vibrations),  such  phenomena  as: 
multiples  of  rotor  frequencies 

inlet  distortions  in  temperature,  velocity  and  flow  direction 
stator-rotor  interactions 
may  be  classified. 

The  second  group  (self-excited  blade  vibrations  or  blade  flutter)  considers 
instead  the  interaction  between  a  blade  movement  and  the  therefrom 
induced  time  dependent  aerodynamic  forces. 

In  a  turbomachine,  flow  phenomena  from  both  these  groups  interfere  with 
each  other  and  create  complicated  unsteady  flow  pattern  tnrougn  tne 
bladings.  This  pattern  can  however  not  be  predicted  with  present 
theoretical  methods,  wherefore  several  idealizations  of  the  time 
dependent  flow  are  normally  made. 


in  the  field  of  self-excited  blade  vibrations,  different  flutter  domains, 
dependent  upon  the  flow  conditions  in  the  machine,  have  been  identified 
throughout  the  years,  both  in  compressors  and  in  turbines.  Due  to  the  nigh 
pretentions- on  modern  jet  engines,  most  of  the  flutter  problems  reported 
have  occured  in  fan  stages  of  these  machines.  Flutter  phenomena  have 
recently  also  become  of  practical  interest  in  transonic  steam  and  gas 
turbines,  especially  at  high  back  pressures.  In  these  operating  conditions, 
the  flow  in  the  blade  passage  is  presumably  transonic  and  eventually 
partially  or  fully  stalled. 

Some  of  the  possible  types  of  blade  flutter  are  represented  schematically 
in  Figure  l.l,  where  a  typical  modern  compressor  chart  is  shown. 

Depending  upon  the  flow  conditions,  the  different  flutter  types  in  the 
Figure  1.1  may  be  characterized  as  (see  for  example  / 2S/-/27/): 

1.  Subsonic/transonic  positive  incidence  stall  flutter  where  the 
compressor  operates  near  the  surge  line  at  either  part  speed  or  near  the 
design  speed.  Although  the  term  “stalled"  is  generally  used,  it  is  not 
settled  that  the  flow  in  this  flutter  region  is  always  separated  over  the 
whole  blade,  or  if  instead  the  flutter  is  associated  with  high  loading. 

2.  Negative  Incidence  transonic  choke  flutter,  at  part  speed. 
This  type  of  flutter  may  appear  when  the  flow  in  a  compressor 
accelerates  through  sonic  transition.  This  flutter  region  is  probably 
associated  with  local  blade  separations  and  unsteady  shock  waves.  An 
important  question  here  is  under  which  circumstances  large  shock 
oscillations  occur. 

3.  Supersonic  positive  incidence  stall  flutter.  The  flow  is 
supersonic  at  the  outer  portion  of  the  blade  and  the  stage  operates  near 
the  surge  limit,  it  is  expected  that  the  flow  has  strong  m-passage 
shocks. 

4.  Unstalled  supersonic  flutter,  at  full  speed  and  with  attached 
flow.  It  can  occur  at  design  point  and  at  higher  or  lower  pressure  ratios 
and  may  limit  high  speed  operations. 

5.  Supersonic  flutter  (type  "A 1 00").  This  type  of  blade  flutter 
has  been  found  to  appear  suddenly  above  a  certain  pressure  ratio,  beiow 
which  no  flutter  is  encountered,  increased  loading  may  shift  this  critical 
pressure  ratio. 

All  of  the  above  mentioned  flutter  limits  are  of  large  practical  interest, 
and  several  investigations,  both  experimentally  and  theoretically,  have 
been  conducted  in  each  region.  It  should  also  be  pointed  out  that  most  of 


the  regions  In  Fig.  l.l  seem  to  have  some  relation  with  the  transonic  flow 
region,  i.e.  they  are  liable  to  be  associated  with  strong  unsteady  shock 
waves.  If  this  unsteady  shocks  have  an  influence  on  the  stability  of  a 
certatn  blading  is  still  a  very  open  question  (the  presented  study  should 
be  wleved  in  this  context  as  far  as  its  shortterm  goal  is  concerned).  An 
unsteady  shock  moving  around  an  a  blade  might  eventually  introduce  some 
phase  lag  between  the  blade  movement  and  the  resulting  pressure 
response,  as  well  as  an  important  change  in  the  amplitude  of  the  lift 
distribution  on  the  blade. 

Presently,  fully  theoretical  models  exist  for  the  prediction  of  unstalled 
self-excited  blade  vibrations  (domains  6,  7  and  8  in  Fig.  1.1),  but  most 
prediction  methods  for  partially  or  fully  stalled  flow  are  heavily  coupled 
with  extensive  experimental  data  sets  to  yield  semi-empirical  models. 
These  flutter  predictions  are  thus  normally  of  proprietary  nature. 

Many  of  the  existing  theoretical  prediction  models  for  unstalled  flow  look 
promising  for  the  computation  of  aeroelastic  forces  acting  upon  vibrating 
blades.  However,  as  the  main  interest  from  the  industry  is  to  compute  the 
unsteady  blade  forces  as  fast  and  cheap  as  possible  in  order  to  avoid 
blade  failures,  only  a  very  limited  number  of  methods  consider  also  the 
fully  unsteady  flow  between  the  blades  although  the  apprehension  of  how 
disturbances  propagate  from  one  blade  to  another  is  essential  to  the 
fundamental  understanding  of  aeroelastic  phenomena  in  blade  rows. 

LZ  Introduction  to  the  Present  Numerical  Methodology. 

It  is  the  purpose  of  the  present  work  to  apply  and  test  a  fully  unsteady 
numerical  method  for,  in  the  future,  predicting  both  the  aeroelastic  blade 
forces  and  the  unsteady  flow  through  vibrating  cascades,  and  thus  to 
contribute  towards  the  comprehension  of  unsteady  physical  flow 
phenomena  in  turbomachines. 

Recently  a  few  studies,  based  upon  the  fully  unsteady  Euler  equations, 
have  been  developed  for  predicting  both  the  aeroelastic  blade  forces  and 
the  unsteady,  unstalled  flow  through  vibrating  two  dimensional  cascades 
/7, 1 5- 1 7/ .  The  study  in  111  snowed  the  feasibility  of  the  approach  for 
two  dimensional  subsonic  flow  through  a  cascade  of  compressor  blades 
with  thin  leading  edges.  However,  the  largest  interest  of  aeroelasticity 
lies  in  the  transonic  flow  region  wherefore  it  Is  considered,  as  an 
extentlon  of  the  work  in  /7/,  to  develop  a  method  with  the  capability  to 


capture  normal  unsteady  shock  waves  (and  later  to  treat  rounded  leading 
edges)  while  maintainig  the  Euler  equations  and  taking  into  consideration 
the  significant  development  of  new  numerical  techniques  during  the  past 
years.  These  techniques  have  mainly  been  directed  towards  the 
consideration  of  the  wave  like  (or  hyperbolic)  nature  of  the  Euler 
equations,  and  performs  the  spatial  differencing  of  the  fluxes  in  an 
“upwind"  (in  regards  to  the  wave  propagation)  manner  (see  for  a  few 
examples  /II -14, 18-23/).  As  pointed  out  in  /23/  these  methods, 
sometimes  denoted  “flux  vector  splitting"  (ref.  /1 1-14/)  and  “flux 
difference  splitting"  (ref.  /1 8-23/),  can  all  be  regarded  as  similar  in 
some  respects  although  the  way  the  fluxes  are  split  are  based  on 
somewhat  more  physical  reasoning  in  /  1 8-23/  and  more  mathematical  in 
/l  1-14/. 

While  this  modern  numerical  methods  seem  to  be  robust  and  reliable  (at 
least  for  steady-state  flows),  a  drawback  which  should  not  be  forgotten 
is  their  large  computational  cost  because  of  the  increased  number  of 
arithmetic  operations  that  are  needed  to  be  solved  per  iteration  step.  The 
methods  are,  and  should  remain,  only  one  aspect  of  a  continuous  ongoing 
effort  to  understand  the  flow  in  turbomachines. 

The  present  study  is  based  on  a  "flux  vector  splitting"  method,  with  the 
application  towards  the  unsteady  flow  through  an  oscillating  cascade. 
Here,  the  method  is  briefly  presented  whereafter  comparisons  with 
existing  numerical  and  exact  (in  one  dimension)  steady  state  solutions  for 
one  and  two  dimensional  channel  flows  are  made.  Test  cases  (in  one 
dimension)  and  numerical  examples  for  unsteady  transonic  flows  with 
oscillating  back  pressure  are  thereafter  given,  as  regards  to  the  methods 
capability  to  accurately  capture  unsteady  normal  shock  waves  in  a  non- 
moving  reference  plane  (i.e.  pressure  oscillations  are  introduced  from 
upstream  or  downstream  and  the  nozzle  walls  and  the  blades  do  not 
move).  A  small  attempt  to  evaluate  how  well  the  numerical  method 
represents  the  expected  physical  phenomena  and  the  implications  of  the 
results  towards  the  blade  flutter  problem  is  given. 


Figures  for  section  l. 


Pres  lure  ratio 


Most  flow 


Fields  6,7  and  8:  Domains  of  application  of  present  fully 
theoretical  subsonic  (6),  transonic  (7)  and  supersonic  (8) 
prediction  models. 


Fig.  1.1:  Different  flutter  types  In  axial  compressors 
(see  for  example  /25/) 


2. 

Nomenclature 

Latin; 

- 

A: 

area 

[m2,-] 

- 

a: 

speed  of  sound 

[m/s,-] 

- 

CP: 

specific  heat  at  constant  pressure 

[m2/(Ks2)] 

- 

Cv: 

specific  heat  at  constant  volume 

[m2/(K$2)] 

- 

e: 

internal  energy/umt  mass 

[m2/s2  -j 

e* 

total  energy/unit  mass= 
=e+q2/2=p/l(*-l)p]+q2/2: 

[m2/s2-] 

®r=iV" 

reference  value  for  energy  and 
entalphy=qr2 

Im2/s2] 

- 

ex-ey: 

unity  vectors 

H 

- 

f: 

frequency 

[HZ] 

- 

F: 

flux  vector  in  £  direction 

[~j 

- 

G: 

flux  vector  in  x\  direction 

f.l 

L  j 

- 

n=e+p/p: 

entalphy/unit  mass 

[m2/S2(-j 

fV 

total  entalphy/unit  mass= 
=h+q2/2=ec+p/p*fyp]/[(-)(- 1  )p]-q2/2 

frrws2,-j 

K: 

reduced  frequency,  based  on  half  chord 
k={w1r}/{2qr}  [-] 

(Note  that  we  here  will  use  q=qr,  i.e. 
the  velocity  used  for  computing  the 
reauced  frequency  is  not  the  inlet,  nor 
the  outlet  velocity.) 

- 

k: 

time  step  indicator 

[-] 

V: 

reference  value  for  length, 
lr=xr=yr=c(=chord) 

[m] 

• 

D: 

pressure 

tN/rn-.-] 

Pamp: 

amplitude  of  sinusoidal  pressure 
fluctuation 

[M/m2,-] 

- 

Pr; 

reference  value  for  pressure 

[N/m2] 

- 

q: 

velocity  vector=ue*+vey 

[m/s,-] 

Qr: 

reference  value  for  veiocify=(RTr)°5 
(Mote  that  qr=speed  of  sound/  v^.) 

[m'S] 

- 

R: 

gas  constant  (*286.7  for  air) 

f  iv»  7  / 1/ 1 

1 1  *  *•*'  W  / 
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V\ 

.M 


% 

- 

R*: 

Riemann  invariant  in  positive 

«!* 

*  i 

_ 

R-:  ' 

coordinate  direction 

Riemann  invariant  in  negative 

H 

coordinate  direction 

H 

3.' 

dimensionless  entropy  =  (s-sr)/cv 

[-] 

»• 

■  § 

- 

t: 

time 

is.-] 

■>  * 

_ 

rcyc!e; 

time  for  a  period  of  oscillation 

Is,-] 

- 

to: 

time  when  oscillation  starts 

Is, -3 

it; 

y 

— 

to: 

old  time  before  next  iteration  step 

ls,-3 

■i 

tr: 

reference  value  for  time=lr/qr 

Is] 

tsh; 

time  for  pressure  disturbances  to 

travel  from  tne  exit  to  the  snock 

is;-3 

~ 

T: 

temperature 

(M 

#v 

Tr: 

reference  value  for  temperature 

rs 

_ 

u: 

velocity  component  in  x  direction 

[m/s,-] 

V 

- 

us; 

shock  velocity  in  x  direction 

[m/s,-] 

$ 

*v 

- 

v: 

velocity  component  in  y  direction 

[m/s,-j 

‘1, 

- 

x: 

Cartesian  coordinate  direction 

[m,-] 

xs: 

shock  position 

[m,-l 

‘5* 

5‘ 

y: 

Cartesian  coordinate  direction 

[m,-] 

$ 

SC3£K‘- 

«,• 

- 

At: 

time  increment  for  one  time  step 

r  1 

IS.-J 

S? 

p: 

density 

[kg/nrp, 

s 

Pr 

reference  value  for  density=pr/(RTr  i 

[kg/m*] 

- 

V 

ratio  of  specific  heats(=i.4  for  air) 

[-] 

£ 

w: 

circular  frequency 

[rad/s] 

, 

l 

curvilinear  coordinate  direction 

r  i 
l"J 

.. 

. 

*< 

_ 

n: 

curvilinear  coordinate  direction 

r  t 

l"J 

l 

«' 

su 

'scripts: 

s 

.« 

1* 

amo: 

amplitude  of  sinusoidal  oscillations 

in  time 

«F 

1  1 

- 

c: 

stagnation  value  in  tne  absolute  frame  of  ref 

1  • 
J 

values  at  downstream  i=ex)t) 

s 

t 


r:  reference  values 

s:  shock 

(>t,  t)A,  ()y  partial  derivative  in  t,  x  and  y  resp. 

infinity  upstream 
+0O:  infinity  downstream 


critical  values 

intermediate  values  in  two  step  integration  method 
dimensionless  values  (only  used  in  ambiguous  context/ 
variables  in  (^q)  plane 


Bold  symbols:  Matrices 


Part  I:  Development,  Implementation  and 
Testing  of  the  Methodology  for  Steady-State 
Flow. 


3.  One  Dimensional  Invlsctd  Steady-State  Flow. 

Before  going  into  the  general  application  of  two  dimensional  flow,  some 
basic  one  dimensional  investigations  are  necessary  to  assess  the  flux 
vector  splitting  method  in  conjunction  with  the  proposed  explicit 
numerical  time  integration  scheme. 

The  one  dimensional  inviscid  gasdynamic  equations  of  motion  can  be 
expressed  in  conservation  law  form  as: 


Wt*Fx=0  (3.0- 1 ) 

where  W  and  F  are  vectors  expressing  the  dependent  variables  and  the 
fluxes: 


w= 

p 

;  F= 

*  « 

pu 

pu 

pu2+p 

Pec 
»  - 

u{pec+p) 

. 

Furthermore,  the  equation  of  state  for  a  perfect  gas  is  given  by: 

ec=p/{('K-1)p}*u2/2  (3.0-3*) 

For  the  present  application,  the  mass,  momentum  and  energy  components 
of  the  flux  vector  F  are  split  into  upwind  parts,  expressing  (in  some 
splittings)  the  "positive"  and  "negative"  eigenvalues  of  F.  Several  methods 
of  performing  such  “flux  vector  splittings"  can  be  found  in  the  literature. 
Out  of  these,  two  have  been  adopted  and  tested  for  the  present 
application.  These  are  the  splittings  according  to  Steqer/ Warming 
(/  l3/,/29/,/l !/)  and  van  Leer  (/28/,/l  1/).  The  basic  idea  behind  any  of 
these  is  that  if  F  can  be  split  into  F=F**F_,  where 

Wt*F*x=0  stable  for  backward  differencing 
Wt*F“x=0  is  stable  for  forward  differencing 

the  same  differencing  can  be  used  for  the  full  equation 
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3.1. 


Flux  Vector  Splitting. 


3.1.1.  Steger/Warming  Splitting. 


in  the  Steger/Warming  approach  the  fluxes  are  split  according  to  the  sign 
of  the  local  eigenvalues  of  the  flux  vector  F.  The  mathematical 
development  of  the  eigenvalues  and  the  splittings  are  given  in  ref. 
/1 3.1 1/  and,  in  a  shorter  version,  in  Appendix  No.  For  briefness,  only  the 
final  results  are  given  here. 

Three  eigenvalues,  X,,  X2  and  X3,  of  the  flux  vector  F  exist  in  one 
dimensional  flow.  These  are  normally  expressed  as: 

X  |  =u 

-  X2=u+a  (3. 1.1-1) 

X3=u-a 


These  eigenvalues  can  be,  depending  on  the  magnitude  of  tne  velocity  u 
(usO,  0<usa,  u>a),  separated  into  positive  and  negative  parts: 

\n«*V*V  (3. 1.1 -2) 

where,  in  tne  original  study  by  Steger/Warming  / 1 3/,  the  X+n,  x_n  are 
defined  as: 

Xn  }/2  (3. 1.1-3) 

Xn  }/2 

These  definitions  where  later  sligthly  modified  by  Steger,  due  to  some 
problems  of  accurately  capturing  the  rarefaction  through  the  sonic  point, 
to  / 29/: 

xV(x„.[\n2.e]°s)/2  (3.1.1  -4) 

X-n=<VlV*e)°5)/2 

where  t  is  a  small  number  in  the  order  of  e=0.04  / 29/.  However,  it  should 
be  pointed  out  that  this  is  a  numerical  trick  (it  introduces  a  numerical 


^+rr^n* 

^"rr^rr 


dissipation  different  from  the  one  created  with  the  use  of  the  relation 
(3. 1.1-3))  only  to  eliminate  a  small  glitch  around  the  some  point,  with  no 
other  apperant  physical  reason  than  creating  continuous  X+,X_  through 
this  point /29, 11/. 

In  whichever  of  the  above  ways  the  X*  and  X*  are  created,  the  flux  vector 
F  can  be  split  into  two  terms  F*  and  F-  / 1 3/,  with: 

F=F+*F"  (3. 1.1 -5) 


where  each  of  the  splitted  flux  vectors  F+  and  F~  can  be  expressed  as 
/ 1 3. 1 1  /: 


F*= 


p/[2*} 


2fy- 1  ]X,  ±u+X2±[u+a]+X3±[u-a] 

hf- 1  jXj  4u2+X2±[u+a]2/2+X3±[u-a]2/2+w 


where 

w={(3-K][X2±+X3i]a2}/{2(if- 1 }}  (3. 1 . 1  -6) 

With  the  definition  (3. 1.1 -3)  of  X*,  it  is  seen  that 

F+=F;  F“=0  if  M=u/a*!  (3,1.1  -7) 

F+=0;  F-=F  if  M=u/as-l 

However,  the  relation  (3. 1.1  -7)  is  not  any  more  valid  if  the  forward  and 
backward  eigenvalues  X±  are  defined  according  to  (3. 1.1-4). 

3.1.2.  van  Leer  Splitting. 

A  somewhat  different  approach  of  splitting  the  flux  vector  F  was  taken  by 
van  Leer  /28, 1 1/.  Here,  the  flux  vector  is  split  into  forward  and  backward 
parts  in  function  of  the  local  Mach  number  8instead  of  in  function  of  the 
eigenvalues^  which  gives: 

For  supersonic  flow  (i.e.  |M|  =  |u/a|*lg  the  relation  (3. 1.1-7)  still  holds: 


(3. 1.1 -7) 
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F+=F;  F~=0  if  M=u/a*1 
F+=0;  F"=F  if  M=u/as-1 


and  for  subsonic  flow,  the  vectors  F+  and  F~  are  expressed  as: 


F*= 


V 

f2* 


where  f1±=tpa{[Mt  1]/2}2 

f2*=fi±-{hH]ui2a}/K 

f3*=f i  ±  {[(^- 1  )u±2ap}/{2(Y2- 1  ]) 


(3. 1.1 -8) 


It  has  been  shown  In  /II/  that  this  splitting  gives  continuously 
differentiable  spilt  fluxes  at  sonic  and  stagnation  points.  This  is  contrary 
to  the  Steger/ warming  splitting  (see  /  1 1/  for  detailed  comparison  of  the 
splitted  fluxes). 

A  recent  review  article  by  van  Leer,  Thomas,  Roe  and  Newsome  /46/ 
indicates  that,  out  of  the  two  splittings  above,  the  van  Leer  splitting  has 
the  smoothest  properties  and  should  capture  a  shock  with  slightly  less 
numerical  grid  points  than  the  Steger/Warming  splitting. 


3.2  Spatial  Differencing 


In  the  flux  vector  splitting  formulation  the  governing  equations  can  he 
expressed, "both  for  the  Steger/ Warming  and  van  Leer  splittings,  as: 

WfFVF-x=0  (3.2-1) 

In  the  original  work  by  Steger/warming  / 1 3/  the  spatial  derivatives  were 
approximated  with  differences  taken  "upwind"  for  both  F*  and  F~,  i.e. 
backward  (in  the  negative  x  direction)  for  the  F+  and  forward  (in  the 
positive  x  direction)  for  the  F~.  In  the  case  of  generalized  coordinates 
this  approach  does  however  not  preserve  "free  stream"  values,  I.e.  it 
introduces  a  dependancy  of  the  transformation  metrics  on  the  resuit  (see 
for  example  /30/).  van  Leer  proposed  instead  (see  for  ex.  / 1 1 /)  to  use  an 
approach,  called  the  "Monotone  Upwind  Schemes  for  Conservation  Laws”  or 
“MUSCL"  approach,  in  which  the  data  is  first  prepared  and  eventually 
limited  before  the  numerical  differences  are  performed.  The  spatial 
derivatives  are  here  approximated  as  centered  differences,  with  values  at 
half  points: 


F±x»AF*j/ Ax={F*  j+ 1  /2_F*i-  \  /2>/Ax 


(3.2-2) 


Furthermore,  the  F±j+j/2  and  F*|-t/2  are  evaluated  with  two  point 
"upwind"  (i.e.  backward  for  F+  and  forward  for  F*)  extrapolation  / 1  i / : 

f 

-|FVi/2=F+(W-j+|/2)  (3.2-3) 

|F-M/2=F-(WV,/2> 

l 

* 

•  W-|»|/2=W|»<t>-| {Wi-W(_|}/2  (3.2-4) 

WV  i  /2=W„  i  i  (W,.  i  -Wi,2>/2 
l 

in  these  expressions,  the  scalar  function  <t>  is  denoted  the  flux  limiter,  it 
has  for  function  to  switch  from  first  order  spatial  accuracy  to  second 
order  at  appropriate  points.  The  value  <t>=l  gives  second  order  spatial 
accuracy  while  $=0  reduces  the  accuracy  to  first  order. 


in  has  been  demonstrated  in  ref.  /  1 1/  (which  is  an  implicit  computational 
code)  that,  due  to  the  fact  that  both  the  AF+  and  AF"  are  evaluated  at  the 
same  compytational  locations  in  the  MUSCL  approach1,  a  method  with  this 
approach  (independent  of  wether  the  Steger/ warming  or  the  van  Leer 
version  is  used  to  split  the  fluxes)  gives  smoother  results  around  a  sonic 
rarefaction  point  than  the  original  Steger/ Warming  spatial  differencing.  It 
also  reduces,  at  ieast  in  ref.  / 1 1 /,  the  overshoots  just  upstream  of  the 
snocK  region  for  one  dimensional  transonic  flows. 


1  Contrary  to  the  original  Steger/ Warming  approach  where  AF+  and 
AF"  are  evaluated  as 

AF*=A°FMA°ackwardF*]={FVF+j- 1 } 

AF-=AfF-=[Af°rwardF-]={F'u ,  -F~i) 
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3.5.  Boundary  Conditions  at  Inlet  and  Outlet 

in  the  present  application  two  different  options  for  treating  the  ooundary 
conditions-at  inlet  and  outlet  were  implemented.  In  the  first  they  are 
treated  in  essentially  the  same  way  as  in  111  where  they  are  implemented 
with  the  "postcorrection  method"  after  Moretti/deNef f  /6/.  In  this 
technique  the  numerical  treatment  of  computational  points  located  on  the 
boundaries  Is  conceived  in  two  steps.  First  the  computation  is  done  using, 
for  the  approximation  of  derivatives  at  the  boundaries,  values  from 
mteneor  points  (and  eventually  extrapolated  therefrom)  only,  in  this 
preliminary  computation  the  boundary  condition  has  not  yet  been  imposed. 
Therefore  the  flow  variables  so  calculated  are  not  correct.  However, 
considering  the  hyperbolic  or  wave-like  nature  of  the  Euler  equations,  the 
Riemann  variables  (a  specific  one  dimensional  combination  of  speed  of 
sound  and  velocity,  see  Appendix  A4)  carried  on  characteristics  impinging 
on  the  boundary  are  expected  to  be  correct.  The  values  of  the  primitive 
variables  (velocity  and  speed  of  sound)  computed  so  far  have  thus  to  be 
changed  by  imposing  the  required  boundary  condition,  while  preserving  the 
value  of  the  Riemann  variable  Just  computed  111 . 

This  updating  is  performed  as  the  final  logic  at  every  time  step. 

In  section  3.3.1  this  postcorrection  technique  is  discussed  for  the  inlet 
and  outlet  boundaries  for  a  one  dimensional  flow, 
in  section  3.3.2  a  second  method,  employed  in  most  calculations  presented 
later,  is  given. 

3.3.1.  First  Met  nod 

The  first  method  for  implementing  the  boundary  condition  is  based  upon 
the  method  used  for  a  solution  with  the  MacCormack  predictor-corrector 
solver  In  /24/,  and  requires  an  extrapolation  of  the  flow  variables  at  the 
boundaries. 

3 . 3. 1  .a .  Extrapolation  of  values  at  the  boundary. 

For  completness,  it  is  here  noted  that  the  points  outside  of  the 
computational  domain  are  extrapolated  with  a  three  point  parabola  from 
the  inner  of  the  flowfield  according  to  the  formula 


where  a,  b  and  c  are  constants.  This  gives  thus,  with  constant  spacing  in  x 
direction: 

If  the  point  (i+1/2)  lies  outside  the  computational  domain: 

a  Ax2  =0.5{f1_2-2f1-1+f1} 

DAx  =0.5-{fj_2-4fj-i+3ft} 
c  =fj 

and  thus 

fM=fi-2-3fi-i+3fi 

U+2~^U-2~%U-\*f>U 

whereafter  the  formula  (3.2-4b)  gives: 

f+i+i/2~^i'*-i'f<^fi*i,^i*i-^i+2^/r2  (3.3. 1  .a-  ic) 

If  the  point  (i-1/2)  lies  outside  the  computational  domain: 

aAx2  =0.5{fj+2-2fj+i+fj} 
bAx  =-0.5-{fi+2-4fi+r3fj} 
c  =fi 

and  thus 

fi-t=fi+2'3fi-H+3fj 
fj-2=3fj+2-8fj+i^6fj 

whereafter  the  formula  (3.4-2a)  gives: 

^"i-1/2=fi-l+<®)~i-1  •{ft-]"ft-2^2  (3.3.1  .a-2c) 


(3.3.  l  .a-2a) 
(3.3. 1  .a- 2b) 


(3.3.  i  .a- la) 
(3.3.1. a- lb) 


3.3  l  b. 


Inlet  Boundary. 


At  the  permeable  inlet  surface  (AB  in  Fig.3.3-l),  the  space  derivatives  m 
x  direction  are  defined  with  an  extrapolation  according  to  the  spar  *'! 
differencing  in  section  3. 3.1. a.  After  the  computation  of  the  whole  flow 
field  (including  the  boundary  points)  the  boundary  condition  is  imposed 
according  to  the  post  correction  technique.  Presently,  two  methods  are 
implemented,  corresponding  to  the  "radiative*  and  "capacity"  conditions 
respectively  ill . 

in  the  first  ("capacity")  technique,  the  flow  is  considered  to  be 
discharged,  from  an  infinitely  large  capacity,  through  a  set  of  small 
nozzles.  In  this  capacity  the  stagnation  temperature  and  stagnation 
pressure  are  constant.  This  leads  to  the  following  two  relationships  for 
the  entropy  and  entalphy: 

s=ln{pc,/pc1r}  =ln{pctO-r).Tclr} 

(=0  if  TC]=PC]=1 )  (3.3. 1  b-  0 

^cl^ci^Y- 1  )=ai  ,new^{(Y~1  ^+^i  ,new2/2 

where  subscript  "new"  indicates  that  the  values  Include  the  correction  for 
the  boundary  conditions. 

in  addition,  a  compatibility  equation  carries  the  information  on  the  left- 
running1  characteristic  from  inside  the  flow  field  onto  the  boundary.  This 
can  be  expressed  as: 

R_ii0id=f(inner  of  flow  field)=known=2aj .oUj/tV-D-qi .old 

(3.3.1.b-2) 

where  subscript  "old"  indicates  that  the  values  do  not  include  the 
correction  for  the  boundary  conditions. 

As  the  value  of  R"  is  assumed  to  be  correct,  we  have 

1  ,old=^~  1  ,new  1  (3.3. 1  .b-3) 


1  Left-running  characteristic  indicates  in  the  present  stud/  a 
characteristic  in  the  negative  coordinate  directions.  Similarly,  the  right- 
running  characteristic  is  in  the  positive  coordinate  directions. 


where 


R”l  ,new=23j  .new^lM  '“"l  .new  (3.3. !  b-4) 

From  the  entalphy  and  compatibility  equations,  aiinew  and  qiinew  can  be 
calculated  as: 

a  i  ,new=iR‘  i  ±(R“  i 2(  1  -K)/2*hc ,  (f  1  )]0-5}(r  i  )/(?♦  |  > 


,new=2ai  new/^K  v3.3.l.b-5) 

whereafter  the  pressure,  density  and  energy  can  be  determined  from  the 
entropy  equation  above. 

The  second  ("non-reflective"  or  "radiative")  technique  to  treat  the 
boundary  is  an  attempt  to  decrease  the  reflections  at  the  inlet. 

The  basic  idea  behind  this  boundary  condition  is  that  in  a  cascade,  in 
contrast  to  an  isolated  airfoil  where  disturbances  propagated  away  from 
the  airfoil  diminish  In  strength  with  the  radial  distance  from  the  profile, 
tne  flow  is  bounded.  The  disturbances  in  an  inviscid  flow  will  thus  not 
always  decrease  in  strength  (for  vibrating  blades,  see  discussions  on  sub- 
and  superresonant  blade  vibrations  /8-10/).  However,  if  the  flow 
upstream  of  the  cascade  is  undisturbed  from  the  left,  the  left-running 
disturbances  will  be  simple  wave  fronts.  If,  after  a  certain  distance, 
these  wave  fronts  become  straight,  it  is  expected  that  the  reflections 
from  an  upstream  computational  boundary  may  be  reduced.  This  way  of 
treating  the  boundary  can  be  visualized  as  an  infinitely  long  duct  with 
undisturbed  flow  upstream  of  the  cascade,  and  it  is  implemented  in  tne 
present  study  in  the  following  way: 

As  in  the  first  technique,  the  compatibility  equation  on  the  left  running 
characteristic  is  considered  to  carry  the  information  from  inside  the  flow 
field  onto  the  boundary: 


R'l  i0id=f( inner  of  flow  field)=known=2ai  i0id/(T“  *  Ml  ,o!d 


(3.3.1.b-6) 

Contrary  to  the  first  technique,  however,  the  boundary  conditions  are  now 
not  imposed  at  the  inlet.  It  is  instead  considered  that  a  right  running 
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characteristic  (coming  from  Infinity  upstream)  carries  information  in  an 
isentropic  way  from  infinity  upstream  to  the  inlet  and  into  the  flow  field. 
This  compatiDility  equation  can  pe  written  as: 

R*  i  .new=^upstream  infinity)=imposed=2a-oor,ew/(K- 1  hq-^new 

(3.3.1. P-7) 

as  the  value  of  R-  is  assumed  to  be  correct,  we  have  as  before 

R~  i  ,oid=R~  1  .new  =23j  .new^Y-'  )“9l  .new  (3.3.  l.b-8) 

The  boundary  conditions  imposed  are  thus  s.*,  and  a.,*,  which  gives 

ai.new={RVR'iHrU/4 

Qi.new=(R>i“R‘|}/2  (3.3.1.b-9) 

it  should  be  pointed  out  that,  as  the  stagnation  values  of  temperature  and 
pressure  are  not  explicitely  imposed  at  the  inlet,  the  "radiative"  boundary 
condition  can  not  always  keep  these  two  variables  at  their  initial  values 
111. 

3 . 3 . 1 .  c .  Outlet  Boundary. 

in  subsonic  flow,  only  one  boundary  condition  may  be  specified  at  the 
outlet.  This  is  usually  taken  to  be  the  static  pressure. 

In  the  present  study,  the  static  pressure  Is  held  constant  (for  steady 
state  flows,  sinusoidal  variation  for  unsteady  flows)  at  the  outlet,  and 
this  boundary  condition  Is  implemented  much  in  the  same  way  as  at  the 
inlet. 

First  It  is  noted  that  the  Riemann  Invariant  commlng  from  the  interior  of 
the  flow  domain  is  calculated  as  above 


R*2,old=2a2,ol<3/(1f_  1  Hz.old  (3.3.1  -C-l ) 

Thereafter  this  value  is  considered  to  be  correct  also  including  the 
boundary  conditions,  wherefore 


Furthermore,  if  the  entropy  calculated  without  considering  the  boundary 
conditions  is  also  considered  to  be  correct  (i.e.  S2,new=S2.oio^  it  is 
possible  to  determine,  from  the  specified  value  P2,new*  the  speed  of  sound 
a2,new  The  velocity  q2,new  is  then  immediately  determined  as 

^2.new=^+2,new_232inew/^K" 1  )  (3.3. 1  .c-3) 

3.3.2  Second  Method 

The  above  method  of  extrapolating  the  flow  variables  to  half  a  grid  point 
outside  the  computational  domain  may  eventually  lead  to  some  problems 
and/or  eventually  introduce  some  unwanted  reflections  at  the  boundaries. 
Another  method  for  treating  the  inlet  and  outlet  boundary  conditions  was 
for  these  reasons  also  developed.  This  second  method  considers  only 
information  from  the  inside  of  the  flow  field  and  combines  the  boundary 
conditions  as  in  the  first  method  (i.e.  Riemann  Invariant  from  the  ieft  or 
constant  stagnation  pressure  and  temperature  at  the  inlet;  given  static 
pressure  at  the  outlet)  with  the  information  carried  along  the  other 
family  of  characteristics  (I.e.  left-running  Riemann  invariant  at  the  inlet; 
right-running  Riemann  invariant  at  the  outlet). 

However,  also  in  this  method  some  special  treatment  is  necessary,  for  the 
2nd  order  spatial  accuracy,  in  the  point  next  to  the  boundary.  In  the 
present  application  this  Is  performed  by  reducing  the  accuracy  of  F+  at 
point  i=2  and  F~  at  point  i=ic-l  to  first  order  spatial  accuracy.  This 
introduces  some  numerical  inaccuracies  at  the  boundaries,  but  if  the  flow 
gradients  are  not  large  (as  at  the  inlet  and  outlet)  or  if  the  mesh  size  is 
reduced  close  to  the  boundaries  (as  close  to  a  leading  edge  of  a  blade),  the 
overall  result  will  still  be  good. 

3 . 3 . 2 .  a .  Outlet  Boundary. 

* 

To  the  Oth  order,  the  following  procedure  can  be  used,  with  reference  to 
Fig.  3.3-2.  The  characteristic  leaving  point  “4"  at  time  t=t0  and  which 
impinges  in  point  "3"  at  time  t=t0+At  is  replaced,  as  in  the  first  method 
above,  by  the  boundary  condition  that  the  static  pressure  is  given,  i.e. 
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p=pe.  Thereafter  the  entropy  is  invariant  along  the  characteristic  between 
points  "2"  and  “3".  Furthermore,  the  right-running  Riemann  invariant  R+  is 
constant  along  the  characteristic  between  points  “1"  and  "3".  The 
following  information  is  thus  obtained: 

P3  =Pe 


R+3  =R*, 

s3  =s2 

wherefore 

a3  =^.p3{(y-D/r)e{s3/T} 
u3  =R+,-2a3/iY-i} 

The  problem  which  then  remains  to  be  solved  is  to  find  the  location  of  the 
points  “1“  and  "2“.  To  the  Oth  order  the  entropy  can  be  obtained  by  a 
simple  extrapolation  from  the  last  point  inside  the  boundary,  i.e. 

s3  =sic- ! 

The  point  "1“  can,  also  to  the  Oth  order,  be  found  by  considering  linear 
interpolation 

{u,*a,)At={xlc-X|} 

{[ui  *ai  Hui_  i  ♦a*- 1  ]}={(ui*ajHui_  i  ♦ai_  1 1Mx! -xj_ }  )/{xpXi_  i } 

where  the  interval  [i- 1  ,i]  is  found  by  considering  when 
{IX|C-xt]/[uraj]-M} 

changes  sign.  From  the  above  relationships,  the  value  of  x,  can  be 
determined  as: 

X|=-((uj-|*at-i]At-xlc-[(Uj*3jHU|-i*3j-i)Jxi-iAt}/ 

{ I  *[(Uj*3|)-(U|-  )]/[X|-X|_  I]} 


wherafter  the  values  of  u  and  a  can  be  found  in  point  i  with  linear 
interpolation  as: 

u,  "={U|-Uj_  j  }{x,  -xj_  |  }/{xrxj_  j  }-Uj_  | 
a,  ={a,-a)_|}{x]-x,-,}/{xrx1.,J+a,-i 

Obviously,  here  a  more  sophisticated  Interpolation  in  order  to  consider  the 
curvature  of  tne  characteristics,  perhaps  aiong  the  lines  as  presentee  in 
ref.  /47/,  would  improve  the  time  accuracy  of  the  results. 

3 . 3 . 2 .  b  Inlet  Boundary. 

At  the  inlet  boundary,  the  boundary  conditions  should  instead  replace  the 
Riemann  invariants  along  the  lines  cr->"3“)  and  C2"->"3“)  (see  Fig.  3.3- 
2),  exactly  as  in  the  first  method  in  section  3. 3.1. a.  Thereafter,  the  point 
"4"  and  tne  values  of  tne  flow  variables  in  this  point  can  De  founo 
according  to  the  treatment  in  section  3.3.2.a  above.  Thus 

x4={(at-rUj-i)At-[(Uj-aj)-(Uj-|-aM)]  Xj-|  At/[xrx|_ibxj=|}/ 

{ 1  +[(urajMUj_I-a1-i)]-At/[xrXj_]}} 

wnereafter 

2a^/(^~  ?  )-u4  =R“«  =R~3 

2a,/('Hbu,  =R*,  =R+3 

and 

u3=0.5{R+i-R-4 

a3=0.5{RVRV/((T-l)/2} 

As  at  the  outlet,  a  more  sophisticated  interpolation  would  not  be  amiss. 


The  quasi  one  dimensional  Euler  equations  can  De  expressed  in  flux' vector 
splitting  form  as  /  I !/ 

W,*(F+*F-)y=-H 


wnere  the  H-vector  takes  into  account  the  variation  of  the  area  A  and  can 
be  expressed  as 


H=(F-P)AX/A 


,  with  P= 


0 

P 

0 


The  numerical  solution  to  this  equation  is  obtained  with  an  explicit 
method,  considering  the  3rea  change  as  a  source  term  with  ak  evaluated 
as  a  centered  difference. 

With  the  above  mentioned  "MUSCL”  approach  / 11/.  an  explicit  second  order 
numerical  scheme  can  be  conceived  as  a  two  step  method  for  the  left  hand 
side  of  tne  above  equation  (see  Appendix  A5  as  well  as  ref  / 12/).  Tnis 
method  is  modified  from  the  original  Steger-Warming  explicit  two  step 
method  for  flux  splitting  / ] 3/  and  C3n  be  expressed  as  / 1 2/: 

First  step:  WVW^-AUACF^+ACFfW/Ax 

Second  step:  Wj^,=0.5-{W*j*W1k-At-(A(F*jM*^fF*ri]/Av> 

where  superscript  indicates  values  at  an  intermediate  time  level  and 
the  spatial  differences  [A(Fj*)'  A(Ft")]  are  expressed  as  centered  between 
points  (i+l/2)  and  (i-l/2): 

A(Fj+)={F*(Wj+  i/2")-F+(Wj-,/2-)} 


A(Fr)={F"(  W.*  i  /2*)-F~(  Wi- 1  /2*)} 


with  the  values  or  the  vector  w  extrapolated  as  follows  /l  1-12/ : 

Wf 1 /2*=Wj+ 1  1  {Wj+ 1  -yH\+2)/2 

Wj+1/2-=Wi>0-r{WrWj-1}/2 

Here  the  scalar  function  $  is  the  flux  limiter  /  1 1/  used  to  switch  from 
first  order  spatial  accuracy  ($=0)  to  second  order  (<t>=l). 
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3-5-  lest  problems  and  Computational  Results 
3.5.1.  First  Test  Problem 

The  first  test  problem  is  the  one  of  a  step  shock  in  a  duct  of  constant 
^ross  section.  Unless  the  exact  solution  is  Initially  specified,  such  a  case 
could  not  be  computed  as  stable.  However,  by  initially  specifying  the 
shock  in  a  certain  position  it  is  possible  to  determine  the  important 
information  of  how  much  the  shock  is  smeared  out  through  the  numerical 
dissipation  innerent  in  the  flux  vector  splitting  method  and  whether  the 
method  can  keep  the  shock  in  the  correct  position,  in  cases  with  no 
influence  of  an  area  change.  Such  a  computation  is  represented  in  Fig. 
3.5. 1-1  for  an  upstream  Mach  number  of  1.2909,  and  for  a  case  where  the 
shock  is  located  exactly  at  a  grid  point  (Steger-Warming  splitting,  with 
61  grid  points).  It  is  seen  that  3  grid  points  lies  in  the  shock  region  and 
that  a  slight  undershoot  is  present  downstream  of  the  shock 
(corresponding  to  0.7 %  of  the  upstream  Mach  number). 

\i  snouia  oe  pointed  out  that,  in  the  presented  calculation,  the  flux  limiter 
function  (0)  was  1,  corresponding  to  a  second  order  spatial  accuracy,  m 
all  grid  points,  except  for  the  calculation  of  the  AF+sh0C(,+3/2,  where  first 
order  was  used  (<J>~ShOck-i=0)  in  order  not  to  extrapolate  the  AF+shoclc+3/2 
over  the  shock  surface.  The  reason  for  this  special  treatment  is  that  F* 
and  F-  are  discontinous  across  a  shock,  although  their  sum  (total  flux=F) 
is  continuous  (see  definition  of  FV1/2). 

Averaging  the  flux  over  the  shock  by,  in  point  (shock- i/2),  defining 
frshock+i/2=^$hock+fr  shock- 1)'2  obviously  eliminates  the  numerical 
dissipation  of  the  flux  splitting  method  and  reproduces  the  exact  solution. 

3.5.2.  Second  Test  Problem 

The  second  test  case  is  the  quasi  one  dimensional  test  case  used  by 
Steger  /! 4/  and  Anderson,  Thomas,  van  Leer  /l  1/  in  connection  with  their 
implicit  flux  vector  splitting  methods,  implemented  with  finite  volume 
approaches.  The  area  of  the  duct  is  given  by,  in  the  present  formulation, 

A(x)  =  1-0.2  {x-0.5Ml.5-x},  with  -0.5sxs1.5. 


(or,  in  the  notation  of  ref.  / 1 1  / 


A(x)  =  1-0.8xa(1-xa)  withOsxAsl). 

The  inlet  and  outlet  boundary  conditions  are  such  that  the  shock  is 
positioned  at  75%  of  the  nozzle  length  /!!/,  with  the  following 
characteristics  of  the  flow: 

filniet  =0.5533 

^upstream  of  shock  =1-2909 
^outlet  =0.5701 

sdown$tream  of  shock  =0.0077 


The  steady  state  results  presented  below  were  obtained  after  a  large 
number  of  iterations,  with  a  total  CFL  (one  dimensional)  in  the  order  of 
CFLtotal=2000-30001 

Fig.  3.5.2- 1  presents  the  Mach  number  and  entropy  for  the  case  of  second 
order  spatial  accuracy  throughout  the  flow  field  (Steger- Warming 
splitting,  with  61  grid  points).  The  shock  is  covered  oy  two  grid  points 
(note  that  the  shock  is  exactly  at  a  grid  point  for  the  exact  solution),  and 
a  sligth  undershoot  (corresponding  to  (Mn-.axundershoor 
^exactJ ^ upstream. of. shock=_ 2% )  is  noted  downstream. 

The  same  results  are  presented  in  Fig.  3.5.2-2  for  the  case  when  flux 
limiting  is  used  for  AF+  in  one  point  downstream  of  the  shock.  Tne  shock 
is  again  represented  by  two  grid  points,  and  the  undershoot  downstream  of 
it  is  somewhat  smaller  (corresponding  to  {Mmax.  undershoot" 
Mexact)/|v1upstream.of.shock=- 1 X).  Similarly,  the  over-  and  undershoots  in 
entropy  are  reduced. 


'  This  value  should  be  compared  to  the  indication  given  in  /I !/  that 
several  hundred  iterations,  with  a  CFL  for  the  there  used  implicit  finite 
volume  method  between  10  and  20,  was  needed  to  obtain  the  stationary 
solution.  The  total  CFL  time  (=CFLtotaP  for  obtaining  the  steady  state 
solution  is  thus  comparable  for  both  methods. 

2  All  the  one  dimensional  steady  state  test  cases  where  run  with  a 
CFl=0.4  per  iteration  step,  i.e.  a  stability  factor  according  to  eg.  (A8-6= 
of  3TAB=0.4. 


If  the  numerical  dissipation  over  the  shock  (introduced  by  the  splitting  of 
the  flux  vector  into  the  two  parts  F*  and  F~)  is  eliminated  Dy  averaging 
the  flux  vector  in  point  (shock*  1/2)  as  FShock* t /2=<FShock*Fihoci.. , ),  the 
shock  is  sleeper  (covered  by  one  grid  point  only),  as  is  seen  in  Fig.  3.5.2- 
3.  However,  as  expected,  the  undershoot  in  both  Mach  number  and  entropy 
is  much  larger  than  in  the  previous  two  cases  (corresponding  to 

^max-undershoor^sxact^  Upstream. of.shock="5%). 

it  is  also  interesting  to  note  that  the  glitch  in  the  rarefaction  through  tne 
sonic  point  as  originally  pointed  out  by  Steger  /14/  (see  for  example  Fig. 

3. 5. 2- 5,  which  is  copied  from  Fig.  10  in  ref.  /  I  A/)  does  not  appear  in  the 
present  method  /  1 1,14/.  This  is  due  to  the  "MUSCL"  extrapolation  of  the 
flux  splitted  values  to  halfway  between  gridpoints  / 1 1  /,  instead  of 
creating  the  flux  differences  in  an  upwind  and  downwind  manner  as 
originally  proposed  by  Steger/Warmlng  /  1 3/. 

in  general,  results  from  the  present  explicit  method  correspond  weil  with 
implicit  results  as  presented  in  ref.  1 1. 

This  can  for  example  be  seen  by  comparing  Fig.  3. 5. 2-6  (copied  from  Fig. 
5,  ref.  /H/)  and  Figs.  3.5.2- la, 2a.  The  undershoots  in  the  shock  region 
seems  to  be  somewhat  smaller  in  the  present  study  than  in  the  results 
presented  in  ref.  /  1 1/,  Fig.  5. 

Results  presented  so  far  have  been  obtained  by,  as  in  / 1 !/,  overspecifying 
the  boundary  conditions  as  the  exact  values  at  the  outlet.  By  instead 
treating  the  outlet  boundary  as  proposed  in  the  section  Boundary 
Conditions,  Method  1”  above  (while  still  overspecifying  the  inlet  boundary 
conditions),  the  results  are  similar  (compare  Figs.  3. 5. 2-2  and  3. 5. 2-4). 
The  only  noticable  difference  is  that  the  Mach  number  in  the  last  point 
upstream  of  the  shock  has  now  decreased  slightly  from  the  exact  solution 
(*0.5%).  Outside  of  the  shock  region,  no  difference  is  noticable  (Fig. 

3. 5.2- 4). 

It  Is  thus  concluded  that  the  present  method  captures  the  shock  in  the 
correct  position  in  a  quasi  one  dimensional  steady  state  flow,  ano  that 
the  smearing  of  the  shock  is  in  general  limited  to  two  grid  elements 
(=3.3%  of  the  nozzle  length  with  61  grid  points).  Furthermore,  the  under- 
and  overshoots  in  the  flow  variables  are  so  small  as  to  be  of  negligible 
importance  for  the  application  under  consideration. 


3.5.3.  Third  Test  Problem. 


The  third  test  problem  is  considered  here  because  of  its  later  application 
in  two  dimensional  flow.  The  quasi  one  dimensional  nozzle  has  an  area 
change  given  by: 

A(x)=0. 1  (5-x)  with  -0.5sxsl  .5  (3.5.3- 1) 

and  the  flow  conditions,  subsonic  along  the  whole  nozzle,  are  governed  by 
an  outlet  Mach  number  of  0.8864  (see  Fig.  3.5.3- 1). 

Tne  quasi  one  dimensional  calculation  as  presented  above  reproduces  rne 
exact  solution  if  6!  grid  points  are  used  (Fig.  3.5.3- la).  A  three  times 
coarser  grid  (21  grid  points)  gives  a  slight  displacement  of  the  isomacn 
lines  (Fig.  3.5.3- 1b,  3. 5.3-2).  As  previous,  it  is  concluded  that  the 
calculation  with  61  grid  points  reproduces  the  exact  solution  within  the 
required  accuracy,  whereas  the  numerical  dissipation  introduced  with  the 
splitting  tends  to  shift  the  solution  from  the  exact  values  if  only  2i  grid 
points  are  used.  In  this  latter  case,  a  non  negligabie  entropy  is  generated 
(Fig.  3.53-3)  and  the  results  from  the  other  flow  variables  are  also 
shifted. 

3.5.4.  Test  of  "Radiative' Boundary  Condition  At  in/et. 

To  find  out  how  large  the  reflection  of  disturbances  from  the  inner  of  the 
flow  field  can  be  expected  to  be  at  the  inlet  boundary  a  one  dimensional 
test  was  run  with  the  "radiative"  boundary  treatment  according  to  section 
3.3.2.P.  Tne  initial  distribution  corresponds  to  a  step  function  tFig  3.5.4- 
1)  and  the  calculation  was  run  until  steady  state.  The  development  of 
pressure  in  the  nozzle  Is  shown  in  Fig.  3.5.4- 1  and  lines  of  constant 
pressure  in  the  (x,t)  plane  is  represented  In  Fig.  3.S.4-2.  it  is  seen  that  a 
certain  reflection  is  present  at  the  inlet  (notice  the  sligth  turn  of  the 
lines  of  constant  pressure  in  time  as  the  disturbance  wave  approaches  the 
inlet),  but  that  the  irregularities  are  not,  in  this  test  case  at  least, 
transported  more  than  a  few  points  into  the  flow  field.  This  is  seen  by 
comparing  Fig  3.5.4-2a  with  Fig.  3.5.4-20  wnere  tne  same  computation  is 
performed  with  the  inlet  boundary  moved  10  grid  points  further  upstream. 
The  isobars  in  both  figures  show  a  similar  behavior.  For  most  practical 
purposes,  this  reflection  is  judged  to  be  of  a  small  nature.  More  tests  of  a 


similar  nature  should  be  performed  also  in  two  dimensions,  especial lv 
with  application  of  the  unsteady  flow. 


Fig.  3.5. 1-1. 

Stepshock  calculation.  First  quasi  one  dimensional  sample  case.  Mach 
number  distribution  along  the  nozzle.  Boundary  conditions  overspecified  at 
inlet  and  outlet. 
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Fig.  3.5.2-1a. 

Mach  number  distribution  for  the  second  quasi  one  dimensional  sample 
case.  Steger/Warming  splitting,  61  grid  points.  Second  order  spatial 
accuracy  throughout  the  flow  field.  Boundary  conditions  overspecified  ar 
•met  and  outlet. 
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Fig.  3.5.2-lb. 

ntropy  distribution  for  the  second  quasi  one  dimensional  sample  case. 
teger/Warming  splitting.  61  grid  points.  Second  order  spatial  accuracy 
throughout  the  flow  field.  Boundary  conditions  overspecified  at  inlet  and 
outlet. 


Fig.  3.5.2-2a. 

Mach  number  distribution  for  the  second  quasi  one  dimensional  sample 
case.  Steger/ Warming  splitting,  61  grid  points.  Second  order  spatial 
accuracy,  flux  limiter  to  first  order  spatial  accuracy  In  first  point 
downstream  of  shock.  Boundary  conditions  overspecified  at  inlet  and 
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Fig.  3.5.2-2b. 

Entropy  distribution  for  the  second  quasi  one  dimensional  sample  case 
Steger/ Warming  splitting,  61  grid  points.  Second  order  spatial  accuracy, 
flux  limiter  to  first  order  spatial  accuracy  in  first  point  downstream  of 
stock.  Boundary  conditions  overspecified  at  inlet  and  outlet. 


Fig.  3.5.2-3a. 

Mach  number  distribution  for  the  second  quasi  one  dimensional  sample 
case.  Steger/ Warming  splitting,  61  grid  points.  Second  order  spatial 
accuracy,  averaging  of  flux  vector  In  point  (shock*  1/2).  Boundary 
conditions  overspecified  at  inlet  and  outlet. 


Fig.  3.5.2-3D. 

Entropy  distribution  for  the  second  quasi  one  dimensional  sample  case. 
Steger/Warming  splitting,  61  grid  points.  Second  order  spatial  accuracy, 
averaging  of  flux  vector  in  point  (shock*  1/2).  Boundary  conditions 
overspecified  at  Inlet  and  outlet. 
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Fig.  3.5.2-4a. 

riach  number  distribution  for  the  second  quasi  one  dimensional  sample 
case.  Steger/Warming  splitting,  61  grid  points.  Second  order  spatial 
accuracy,  first  order  in  first  point  downstream  of  shock.  Boundary 
condition  overspecified  at  inlet,  constant  pressure  at  outlet. 
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Fig.  3.5.2-^b. 

Entropy  distribution  for  the  second  quasi  one  dimensional  sample  case. 
Steger/Warmlng  splitting,  61  grid  points.  Second  order  spatial  accuracy, 
first  order  in  first  point  downstream  of  shock.  Boundary  condition 
overspecified  at  inlet,  constant  pressure  at  outlet. 
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Fig.  3.5.2-Vc. 

Density  distribution  for  the  second  quasi  one  dimensional  sample  case. 


Steger/Warming  splitting,  61  grid  points.  Second  order  spatial  accuracy, ; 
first  order  in  first  point  downstream  of  shock.  Boundary  condition 
overspecified  at  inlet,  constant  pressure  at  outlet. 
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Fig.  10.  Computational  comparison  of  \ *  formulations. 


Ftg.  3.5.2-5  (=Flg.  10  In  ref.  /I A/). 

Illustration  of  "glitch"  at  sonic  rarefaction  for  quasi  one  dimensional 
flow. 


Distance  along  nozzle 


(b)  Steger-Waratng  splitting  («  •  0) 


Figure  S.  -  Mach  number  distribution  along 
~  transonic  nozzle  using  MUSCL>type 
differencing. 


Fig.  3. 5. 2-6  (=F1g.  5  In  ref.  /!!/). 

Mach  number  distribution  from  ref.  /  I  t/.S 
Please  note  the  different  definitions  of  nozzle  length. 
Ref. /11/:  OsXsl;  xs=0.75 

Present  study:  -0.5sxsl.5  xs=1.00 


Ftg.  3.5.3-1a. 

Mach  number  distribution  for  the  third  quasi  one  dimensional  sample  case 
Steger/Warming  splitting,  61  grid  points.  Second  order  spatial  accuracy 
first  order  in  first  point  downstream  of  shock. 


Fig.  3.5.3-lb. 

Mach  number  distribution  for  the  third  quasi  one  dimensional  sample  case 
Steger/ Warming  splitting,  21  grid  points.  Second  order  spatial  accuracy 
first  order  in  first  point  downstream  of  shock. 
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Fig.  3. 5. 3-3. 

Entropy  distribution  for  the  third  quasi  one  dimensional  sample  case. 
Steger/ Warming  splitting,  21  and  61  grid  points.  Second  order  spatial 
accuracy,  first  order  in  first  point  downstream  of  shock. 
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Fig.  3.5.4- 1  . 

Pressure  along  the  nozzle  at  different  time  steps,  fourth  quasi  one 
dimensional  test  case,  "non-reflecting"  inlet  Ooundary  condition  according 
to  method  l  (section  3.3. 1.b). 
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Fig.  3.5.4-2D. 

isobars  In  the  (x,t)  plane  for  the  fourth  quasi  one  dimensional  test  case. 
Initial  pressure  distribution  according  to  Fig.  3.5.4- 1.  "Non-ref lective" 
Inlet  boundary  condition  according  to  method  I  (section  3.3.1.b).  Same 
calculation  as  In  Fig.  3.5.4-2b  but  computational  inlet  moved  10  grid 
points  upstream. 


fi; 


4.  Two  Dimensional  Steady  State  Invisctd  Flow. 


The  rotational  invlscid,  unsteady,  two  dimensional  flow  (without  mass 
forces,  neat  transfer  and  exchange  of  energy),  considering  a  perfect  gas, 
can  be  described  in  differential,  conservative,  form  by  the  continuity, 
momentum  and  energy  equations  (see  for  example  /3/): 


Continuity: 
Momentum  in  x: 

in  y: 

Energy: 


pfKpu)x*(pv)y 

(pu) t+(pu2+p)x*(puv)y 

(pv) t+(puv)x+(pv2+p)y 
(p[e+q2/2])t+(pu[e+p/p+ 


=  0  (4.0-1) 

=  0  (4.0-2) 

=  0  (4.0-3) 

q2/2])x+(pv(e-*-p/p+q2/2])y=0 


(4.0-4) 


as  well  as  the  equation  of  state: 


Perfect  gas:  p/p=RT  (4.0-5) 

The  energy  equation  can  also  be  written  in  the  form: 

Energy:  (pec)t+(pncu)x+(phcv)y  =  0  (4.0-6) 

where  ec  and  hc  are  the  total  energy  and  entalphy  per  unit  mass 
respectively,  with 

ec  =p/{phHD+q2/2 
hc  =ec+p/p  =rp/{phH]}+q2/2 

These  equations  are  Identical  in  dimensionless  form  if  the  following 
reference  values  are  used  (see  Appendix  A1  for  validation): 

Pr>  *r» 

pr=pr/(RTr) 

xr=yr=c(=chord) 

qr={RTr}° 5  (Note  that  qr=speed  of  sound//?.) 

tp=c/ qr 

er=nr=qr2 

s=(s-sr)/cv 


From  here  on,  all  values  are  thus  expressed  in  dimensionless  form  unless 
stated  otherwise. 

With  these  reference  values,  the  relationship  between  the  flow  variables 
p,  q,  p,  ec  and  a  (=speed  of  sound),  M  (=Mach  number),  T  ^temperature)  is 
given  in  Appendix  A3. 

The  equations  (4.0-1  to  4.0-4)  above  can  be  expressed  in  vectorform  as 


with 


W= 


Gy  = 

=  0 

(4.0-7) 

P 

pu 

pv 

pu 

;  F= 

pu2+p 

;  G= 

puv 

pv 

puv 

pv2+p 

Pec 

pu{ec+p/p> 

pv{ec+p/p} 

J 


J 
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4.1.  Governing  Equations  In  Generalized  Coordinates. 


The  equation  system  (4.0-8)  can  be  mapped  from  the  physical  plane  of 
reference  (x,y,t)  into  a  computational  plane  of  reference  ($,ti,t)  with  a 
general  transformation: 


H(x.y.t) 

n=n(x,y,t) 

T=t 


(4.1-1) 


Viviand  /A/  has  shown  that  eq.  (4.0-8)  can  stay  in  strong  conservative 
form  with  a  transformation  of  the  above  type.  The  eq.  (4.0-8)  is  then 
transformed  Into  (see  Appendix  A2): 


wyFyGyo 


(4.1-2) 


with 


W* 

F* 

G' 

D 


D-<W 

D-%>KxF^yG} 

D-HritW+nxF^y6) 

SxHyMx  =  l/^yn-xny?} 


(4.1-3) 


where  D  is  the  Jacobian  determinant  of  the  transformation  (4.1-1)  and  the 
metric  terms  are  related  to  the  derivatives  of  x,  y  and  t  by  (see  for  ex. 
/31 ,  eq.  (4)/): 


$x 

Sy 


-xnD 

~Xr^x_yT^y 


(4.1 -4a) 


lx 

% 

nt 


L 


-y  tP 

x^D 

-xTnx-yrny 


(4.1  -4b)' 


f  *  --W  ^ M  *  1 "C  *  •f' .  ■'*  1  ■»"« 


~7t)’ 

z 


The  Jacobian  determinant,  In  the  center  of  an  element,  corresponds 
physically  to  the  inverse  of  the  cell  area  (Fig.  4.1-1). 

In  component  form,  eq.  (4.1-3)  can  be  written  as: 

P 


W*=  D-l 


pu 

pv 


(4.1 -5a) 


|  Pec 
l 


F'=  D"1 


^tP^xPU^ypV 

$tpU+$x{pU2+pKypUV 

^tpV^xpUV^y{pv2*p} 

$tpec+$xpu{ec+p/pKypv{ec+p/p} 


(4.1 -5b) 


G‘=  D-1 


ntP+Hxpu+nypv 

i]tpu+iix{pu2>p}+nypuv 

TltPv+TlXPuV+T1y{pV2+p} 

ntec+nxpu{ec+p/p}+nypv{ec+p/p} 


(4.1-50 


in  the  first  part  of  the  present  study  we  treat  only  non-moving  grids, 
wherefore  the  terms  ^  and  nt  are  zero.  The  expressions  (4.1-5)  are  then 
reduced  to  (see  for  example  ref.  /II/,  eq.  (35)): 


W'=  0" 


P 

pu 

pv 

Pec 


(4.1 -6a) 


F‘=  D- 


^XpU^yPV 
^X{pU2>p}^yPUV 
^XpUV^y{pV2  +  p} 


(4. 1  -6b') 


”7 

/ 


K 


^xpu{ec*p/p}+^ypv{ec*p/p} 


or 


F'= 


ynPu-xnpv 

yn{pu2+p}-xnpuv 

ynpuv-xn{pv2+p} 

ynpu{ec+p/p>-xnpv{ec+p/p} 


or 


F'= 


f, 

f‘i^ynp 

f,v-xnp 

f'i{ec+p/p} 


as  well  as 


G'=  D-» 


G‘= 


TlxpU+HypV 

r]x{pu2+p}+riypuv 

r)xpuv+r|y{pv2^p} 

Hxpu{ec+p/p}+i]ypv{ec+p/p} 


-y^pu+x^pv 

-y^{pu2+p}+x^puv 

-y^puvfx^pv2+p} 

-y*pu{ec*p/pbx^pv{ec+p/p} 


G‘= 


g'i 

gvj-y^p 

g’,v*x^p 

g'i{ec+p/p) 


(4.i-6tn 


(4. 1  -6P‘") 


(4.1-60 


(4.1-60) 


(4.1-60) 


D-1 


Area  = 


{x*yn-xny*} 

{[x4-x8][y6-y2Hx6-x23[yry8]}/(U4-^l[n6-n2]} 

-0.25{[x5^x3-xrx1][y5+yry3-y,]-[x5*x7-x3-x1][y5+y3-y7-y1]} 

/^s4“^8K136~n2^ 

0.25-{x5y5+X5y7-X5y3-x5y^x3y5+X3y7-X3y3-x3yrx7y5-x7y7 

♦x7y3+x7y, -x, y5-x, y7+x, y3+x, y , 

-X5y5-X5y3+x5y7+x5yrx7y5-x7y3^x7y7^x7y1+X3y5+x3y3 

■x3y7'X3y1+x1y5*x1y3-x1y7-x1y1}/{[^-^8][T|6-r|2]} 

0.5-{x5[y7-y3]*X3ly5-y1]*x7lyry5]*X|[y3-y7]>/{Ur^8llTi6-^2J} 

0.5{[x5-x,  Ky7-y3Hx3-x7][ys-y,  MlMsKne-y) 

0.5  {Ix5-x, ][y7-y3]*lx3-x7l[y5-y ,  ]}  if  A$=AiH 

[y7-yi  3Cx7-x,  ]/2+[y5-y7][x5-x73/2+[y7-y !  ][x5-x7] 

-fx3-x,  J[y3-y  t  ]/2-fxs-X3]fy5-y3]/2-rxs-X3J[y3-y ,  ] 

0  5-{y7x7-y7x, -y, x7*y, x, +y5x5-y5x7-y7X5+y?x7 

+2y7xs-2y7x7-2y,  x5+2y ,  x7-x3y3*x3y , +x,  y3-x,  y , 

■2x5y3>2x5yj+2x3y3-2x3y1-x5y5+x5y3+x3y5-x3y3} 

0.5  {y7xs-y7x, +y ,  x7-y , X3^y5X3-y5x7+y 3x, -y3x5) 
0.5{y7[xs-x,]-y3(x5-x,3+y5[x3-x7j-y,[x3-x7]} 
0.5*{lx3-x7]ly5-yi  Hx5-x,  ][y7-y3J) 


Fig.  4.  1-1:  Geometrical  interpretation  of  the  Jacobian  determinant 


<4.2  Grid  Construction. 


A  grid  for  numerical  calculations  can  be  constructed  in  several  different 
ways,  for  -example  in  a  finite  difference  or  finite  volume  manner.  The 
authors  previous  work  has  consisted  of  methods  when  the  grid  points  are 
designed  to  coincide  with  the  boundaries  of  the  computational  domain,  in 
the  present  study,  some  -problems  appeared  in  the  implementation  of  such  a 
grid  in  connection  with  the  imposed  boundary  conditions  in  the  flux  vector 
splitting  formulation.  To  this  end  the  numerical  grid  points  where  instead 
constructed  as  in  a  finite  volume  approach,  i.e.  the  grid  points  are  in  the 
cell  center  and  the  cell  interface  coincide  with  the  computational 
boundaries.  As  both  approaches  where  used,  both  will  briefly  be  explained. 


The  grid  choosen  first  was  one  where  the  grid  lines  coincide  with  the 
boundaries  of  the  computational  domain,  In  agreement  with  the  author's 
previous  work  /24/  (see  Fig.  4.2. i-l).  In  this  previous  work,  the  unknown 
quantities  (outside  the  computational  domain)  have  been  parabolic ly 
extrapolated  (according  to  section  3.3.1. a).  This  gave  good  results  in  /24/( 
where  numerical  method  by  MacCormack  /48/  was  employed  together  with 
the  postcorrection  technique  after  Moretti  /6/.  However,  with  this  grid 
construction,  the  present  Implementation  of  the  flux  vector  splitting 
method  showed  large  convergence  problems  for  H-grids. 

The  reason  therefore  is  beleived  to  be  found  at  the  boundaries. 

In  a  grid  as  in  Fig.  4.2.1- 1  the  implementation  of  the  boundary  conditions 
is  not  straightforward  in  the  T)USCL"-approach  of  the  flux  vector  splitting. 
This  is  immediately  seen  from  Fig.  4.2. 1-1,  where  it  is  noted  that,  for 
computation  of  the  flow  variables  in  point  j,  the  flow  variables  in  point 
“J-1/2"  are  needed  (as  (Gjiy^Gj+i^-Gj-^J/Ay).  This  does  not  pose  any 
problem  for  the  flux  G~j_i/2  which  is  cunstructed  with  information  from 
points  “j“  and  "j*r,  i.e.  from  known  values.  However,  G*j_i/2  should  be 
constructed  with  information  from  points  ”j-r  and  "j-2",  according  to  eq. 
(3.2-4).  This  implies  an  extrapolation  of  two  points  outside  the  boundary 
for  all  the  flow  variables. 


The  second  method  to  construct  the  grid  is  instead  to  put  the 

computational  points  in  the  middle  of  a  computational  element  (Fig.  4.2.2- 
l).  In  this  case  the  computation  of  the  flow  variables  on  line  "J"  consits 
of,  as  before,  (Gj]y={Gj+i/2-Gj-i/2)/Ay.  However,  “j-l/2"  coincides  now 
with  the  wall,  wherefore  the  boundary  condition  of  no  mass  flow  through 
the  wall  gives 

»  ^  ^ 

pv  0 

Gj-i/2”  pvu  =  0 

pv2+p  p 

p[ec+p/pj  0 

i  J 

it  is  therefore  not  necessary  to  perform  the  flux  vector  splitting  at  the 
wall.  It  should  however  be  noted  here  that  this  procedure  introduces  a 
different  numerical  dissipation  at  the  wall  than  in  all  other  points.  The 
reason  for  this  is  obviously  that  the  total  flux  in  point  "j-l/2"  is  without 
the  numerical  dissipation  the  flux  vector  splitting  introduces  in  the 
general  points.  However,  if  the  cells  are  small  close  to  the  blade  surfaces, 
the  numerical  dissipation  Is  small  anyway,  wherefore  no  problems  are  to 
be  expected  for  a  grid  that  clusters  points  close  to  the  blade  surface. 

Thus,  instead  of  extrapolating  all  flow  variables  to  two  points  as  in 
section  4.2.1,  it  is  only  necessary  to  determine  the  pressure  1/2  point 
away  from  the  known  value. 

The  value  of  the  pressure  in  point  “J*  can  be  found  by  several  methods,  for 
example  by  extrapolation; 

For  the  lower  boundary 

Oth  order  extrapolation;  Pj-i/2,i  =Pj,i 

1st  '  -  :  Pj- 1/2.1  =3pj/2-pj.|,j/2 

2nd  '  ‘  :  Pj-1/2.1  =(5P),|.2'l8P).|,i*2IPji|!/8 

For  the  upper  boundary 

Oth  order  extrapolation:  pj+ 1/2,1  =pjj 

1st  "  B  :  Pj*1/2,i  =3pjt|/2-pj-!,|/2 


Better  still,  in  the  case  of  an  orthogonal  grid,  is  to  solve  the  momentum 
equation  normal  to  the  wail  and  to  obtain  the  pressure  therefrom  (see  for 
example  ref.  /n/). 

For  the  steady  state  sample  cases  in  section  4.6,  only  the  simple  Oth  order 
extrapolation  was  employed. 


43.  flux  Vector  Splitting 

43.1  Cartesian  Coordinates. 

4.3. 1  .a.  Steger/ Warming  Splitting. 

As  in  one  dimensional  flow,  the  fluxes  in  the  Steger/ Warming  approach 
are  split  according  to  the  sign  of  the  local  eigenvalues  of  the  flux  vectors 
F  and  G. 

The  equation  (40-8)  has  four  eigenvalues,  Xjx,  X2\  X3X  and  X4\  of  the 
flux  vector  F  in  the  x  direction.  Similarly,  four  eigenvalues  of  the  the  flux 
vector  G  exist  in  the  y  direction  (Xj/,  \2y,  X3y  and  X4y).  The  eigenvalues 

are  normally  expressed  as: 

- 

X,*=u 

X2x=u 

■  x3x=u+a  (4.3. 1  .a- 1  a) 

x4x=u-a 


x,y=v 

x2y=v 

-  x3y=v+a  (4.3.  i  .a-  i  p) 

x4y=v-a 

These  eigenvalues  can  be,  depending  on  the  magnitude  of  the  velocities  u 
and  v  (u,vsO,  0<u,vsa,  u,v>a),  separated  into  positive  and  negative  parts  as 
in  one  dimensional  flow: 


Xn=X+n+X  r 


(4.3. 1  a-2) 


where,  in  the  original  study  by  Steger/Warming  /l  3/,  the  X+n,  x*n  are 
defined  as: 


As  in  one  dimensional  flow,  the  flux  vectors  F  and  6  can  be  split  into  two 
terms  F+,6*  and  F~,G~  / 1 3/.  with: 


F=F**F- 

G=G+*G- 


(4.3.1.a-5) 


where  the  splitted  flux  vectors  F+  and  F~  can  be  expressed  as  (see 
Appendix  A7  and,  for  example,  ref./ 1 1/,  eq.  (28)): 

1 


P±  =  7- 
2y 


f2(^1A,«  +  X3x± 

i 

+  X4x± 

2(r1)X,  V*  +  W* 

>  W 

l2(ri)vX,x±  +  vx3<± 

+  vX4x± 

i  X.  x2+v2 

UrllqV^Pz^* 

X4x2+v2 
+  [  2 

whith 


„♦  (S-yHV^-fOa2 

w1  =  - - - 

2(Y-1) 

and  the  splitted  vectors  G*  and  G"  as  (Appendix  A7): 

f2(y1)X," 


(4.3.1  .a-6a) 


,«  +  u!Uy± 


Qt_  P_  ^(T1  )uXi 

2lf  l2(y1)>.,yX/± 


+  uX4y± 


(r-i)q 


+  X/X/*  I 

n2^  y2  u2+x y2  I 

v* + -  wy±  j 


where 


w 


y± 


(3-Y)(^3yt->-^4yt)a2 


(4.3. 1  a-6b) 


2(yi) 

As  for  one  dimensional  flow,  it  is  seen  that 


F+=F;  F_=0  ifMX=u/asi 
F*=0;  F~=F  if  Mx=u/as- 1 


(4.3.l.a-7a) 


G+=G;  G“=0  ifMy=v/ail 


(4.3.1.a-7b) 


e-=0;  G_=G  if  ny=v/ai-l 


However,  as  for  one  dimensional  flow,  the  relations  (4.3. 1-7)  are  not  any 
more  valid  if  the  forward  and  backward  eigenvalues  X*  are  defined 
according  to  (3. 1.1 -4). 


4.3.1 .b. 


van  Leer  Splitting. 


As  in  one  dimensional  flow,  van  Leer  proposes  to  split  the  flux  vectors 
into  forward  and  backward  parts  in  function  of  the  local,  one  dimensional, 
Mach  number  instead  of  in  function  of  the  eigenvalues,  which  gives: 


For  supersonic  flow  (i.e.  |M*|*1  and/or  |MY |x  1 )  the  relations  (4.3.1.a-7) 
still  hold: 


F*=F;  F'=0  if  Mx=u/a*1 
F+=0;  F"=F  ifM*=u/as-l 


(4.3.l.a-7a) 


G+=6;6-=0  ifMy=v/a*1 
G+=0;  G"=G  if  My=v/a*-l 


(4.3.1.a-7b) 


and  for  subsonic  flow,  the  vectors  F+,G*  and  F~,G_  are  expressed  as 
/1 1,28/: 


F*=  f,* 

V 

u* 

»  a 


(4.3.1.a-8a) 


where 


f,*=ipa{[Mxi  l]/2}2 

f2±=V-{hH)ut2a}/)f 

f3±=f|±.v 

V=f ,  *  ■[{[(-*- 1  )ui2a]2}/{2[^2- 1  })>v2/2] 
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G*= 


where 


gi* 

92* 

93* 

9<* 


g]±=±pa{[My±!]/2}2 

g2±=g,±.u 

93±=gi  ^  ]v±2a}/> 

g<*=g,  *•[{!(?-  l  )v*2a]2}/{2(K2- 1  ]}*u2/2] 


(4.3.1  .a-8b) 


4.3.2  Generalized  Coordinates. 


I 


The  original  flux  vector  splitting  study  by  Steger/Warming  /  1 3/  contained 
a  splitting  In  generalized  coordinates  which  was  different  from  the  one  in 
cartesian  coordinates,  and  consisted  of  a  mixture  of  %  and  t\  derivatives  in 
both  the  F*  and  G‘  fluxes  (see  section  4.1  for  definition  of  F‘  and  G‘). 
Anderson/Thomas/van  Leer  /II/  proposed  instead  to  locally  rotate  the 
fluxes  into  a  “body-fitted"  coordinate  system,  first  for  the  F  flux,  while 
treating  the  G*n  as  a  source  term,  and  thereafter  the  G'  flux,  while 
treating  the  F’^  as  a  source  term.  This  procedure  is  explained  in  detail 
below  (see  ref.  / 1 1/). 

The  equation  (4.1-2)  can  be  multiplied,  from  the  left,  by  a  transformation 
matrix  T*: 


(TtWV(T*r)^-T*G' y(T*)TW-(TtyF’ 

which  by  defining 

■W  =T$W* 

F  =TSF‘ 


(4.3.2- 1) 


(4.3. 2-2) 


reduces  to 


('W)r+(F)^=-T^G,rl*(T^)TW,+(T^)jF' 


(43.2-3) 


Splitting  of  (F)^,  yet  undetermined,  gives: 
(W)T+(F*+F-)*=-T*Gy(T*)TW-*(TtyF- 


(43.2-4) 


whereafter  multiplication,  from  the  left,  with  T*-1  (=the  inverse  of  T*) 
gives; 


TS- '  (^'t-T?tW  bit- 1  (F*+F-)p-T*- 1  (T*Gy*T*- 1  (T*)<F‘ 


(43.2-5) 


I  (T*W  V+(T$~ 1  F**T^ 1  F-)r(T*“ 1  )^(F**F-)= 
=-T^“,(T^G,^)+T^",(T^)^F‘ 


(43.2-6) 


I  W‘T*  (T$-  1F++T*-  5F-)r(T^- 1  )*(F++F‘)= 
=-IG*rj+T^~HT^)^F' 


(43.2-7) 


where  I  identity  matrix= 


10  0  1 
0  10  0 
0  0  10 
0  0  0  1 


This  implies  furthermore  that 


IWy(F,+*F-)^IG’n= 


where 


F'*  =lTt)-?T* 

F-  =(T^)_,F' 


(4  3.2-8) 


(43.2-9; 


(4  3.2-10) 


Thereafter,  the  RHS  of  equation  (4.3. 2-9)  becomes,  as  F=T*F’ 


(4.3.2-11) 


T?-1(T^)?F,>(T^1)j(F^T-)=TM(T^)?F^(T^,),(T^F,)= 

=(T^_,T^)^F’=I^F'=0 


Equation  (4.3-9)  becomes  thus: 

WT*(F,t*F,-)^6,J1=0  (4  3.2-12) 

which  is  identical  with  the  original  equation  (4.1-2). 

It  is  interesting  to  note  that  the  term  G*n  has  not  been  affected  by  the 
above  manipulation.  Thus,  it  is  possible  to  perform  a  similar  procedure 
also  for  tms  flux.  Tms  gives: 

wy  (F**+F- V(G,+*G‘-)n=0  (4.3.2- 1 3) 


with 

G'MTUr1^  (4.3.2-14) 

G-*(H)- 113- 


Note  that  no  assumption  has,  up  till  now,  been  made  concerning  how  the 
flux  splitting  should  be  performed. 

if  now  the  rotation  matrices  T*  and  Pi  are  choosen  in  a  special  way,  the 
transformed  fluxes  T  and  15  can  be  of  the  same  form  as  the  original 
equations,  wherefore  any  of  the  splittings  (4.3.i.a-6a)  and  (4.3. i .a-6b),  as 
developed  for  the  cartesian  coordinates,  can  be  applied  directly. 

To  obtain  such  a  form  for  T  and  75,  define  the  rotation  matrix  T$,  in  the 
case  of  a  non-moving  grid,  as: 


T*= 


with 


0 

0 


0  0 

C*|  C*2 

-C*2 

0  0 


0 

0 

0 


(4.3.2-15) 


c?i=V“x2+5y2}0  5=yn'{yn2*xTi2}0  5 


<  4  7  o_  1  c.\ 


nm  w^rwr:  vum 


^njTTTwvr^girgirginnn<^v w\ jntww.v'  vi  *-n nr*  -^«r vww  i»  v»  ^  u’lt  ^ 


C?2=V^x2^y2)°  5=-><n/ W^n2)0'5 


This  rotation  gives  then: 


W=T^W’=D"' 


P_ 

pu^ 

pv* 

Pec 


L  J 


T=T^F*=D-’Ux2^y2}0  5 


{yn2*xn2}° 


.5 


ptf 

pu*  u^p 

pv* 

pu^[ec>p/p] 


ptf 

put  &*o 

pvtut 

pu*[ec+p/p] 


o- 


(4.5.2- !  7a ' 


(4.3.2- 17b) 


where 

uM^u^yv}/{^^y2}0.5={ynu-xnv}/{yn2+Xj12}0.5  (4.3.2-  17C) 

^{-?yu^xv}/Ux2^y2}0.Mxnu*ynv}/{y^xn2}0.5 

are  the  rotated  velocity  components,  ttf  is  the  velocity  normal  to  a  imp  or 
constant  %,  representing  the  scaled  contravariant  velocity  component 
/ 1 1/,  and  v*  is  normal  to  u*  (Fig.  4.3.2-la). 

As  the  eqs.  (4.3.2- 1 7a,  17b)  are  in  the  same  form  as  the  eqs.  (4. 1 -6a, 60) 
for  the  cartesian  coordinates,  splittings  identical  to  the  ones  in  the 
cartesian  system  can  be  used  for  splitting  of  the  flux  vector  F‘  in 
generalized  coordinates. 

in  a  similar  way,  let  us  define  for  the  rotation  matrix  Tn  (non-moving 
grid): 


0  0 

-cn,  cn2 
C^2 

0  0 


with 


(4.3.2-18) 


<*,  _ny//^nx^+ny^^  ^-Xjj/{X;»2+y£2}0 ,5 
0^2=1]%;/  (rjvt+r]y2)0i5=-y^/  (x^2+y^2)0.5 


This  rotation  gives  then: 


■W=Tnw=D-i 


P 

pun 

p\7n 

Pec 

(.  J 


pvH 

pUn  vn 

G=TnG,=D',  {r]x2+r]y2}0 5  pvn  vH+p 

pvn[ec+p/p] 

L  J 

*■ 

pvH 

pijn  vn 

=  {y?2+X^2}0-5  pvH  vH+p 

pvn[ec+p/p] 


/here 


uri-(-nyu+r]xv}/{r]x2+r|y2}0  5={-x^u-y^v}/{x52+yt2}0.5 

vn={r)xu*r)yv}/{nx2+iiy2}0  5={-y?u^v}/(x^2)0  5 


(4.3.2-19) 


(4.3.2-203) 


(4.3.2-20D) 


(4. 3. 2- 20c) 


r'  r* 


rwwwirwnwKM jqnwrwTtnrjyuwvru w~rr*  w  uvw /■« .  ^ 

Ch 
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are  the  rotated  velocity  components,  v1)  is  the  velocity  normal  to  a  line  of 
constant  rj,  representing  the  scaled  contravariant  velocity  component ,  and 
Un  is  normal  to  vn  (Fig.  4.3.2- lb). 

As  the  eqs.  (4.3.2-20a,20b)  are  in  the  same  form  as  the  eqs.  (4.1 -6a, 6c) 
for  the  cartesian  coordinates,  splittings  identical  to  the  ones  in  the 
cartesian  system  can  be  used  for  splitting  of  the  flux  vector  G’  in 
generalized  coordinates. 

After  performing  the  splittings  of  7  and  "5,  these  fluxes  have  to  be 
rotated  back  into  the  F’  and  G‘  fluxes  with  the  inverse  of  the 
transformation  matrices  T*  and  Pi  (see  eqs.  (4.3.2-10,14)).  These  can  be 
expressed  as: 


TH= 


I 

0 

0 

0 


0 


<*1 

C$2 

0 


0  0 

-c*2  0 

c$,  0 

0  I 


(4.3.2-21) 


Tn-i=m= 


10  0  0 

o  -cn,  cn2  o 

o  cv2  c^  o 

ooo  i 


t  A  ~7  O1 


V 


mmsmm 


Definition  ot  us 


b<  Definition  of  u’i,  v1* 


Ftg.  4.3.2- 1. 


Definition  of  contravariant  velocities 


4  a  /.  Boundary  Conations  at  the  Blade  Walls. 

4.4.  i. a.  Grid  Points  Situated  on  the  Walls. 

in  the  C-qnd,  the  treatment  of  the  boundary  conditions  at  the  walls  (no 
mass  through  the  wall)  are  straightforward  as  a  mirror  image.  This  is  the 
case  as  the  grid  can  be  continously  generated  in  one  point  inside  the 
ooundarv  (as  tne  grid  is  then  orthogonal  to  the  boundary)  and  tnus  the 
fluxes  can  be  split  in  the  same  way  as  in  the  interior  points.  Then,  me 
consideration  of  a  mirror  image  immediately  gives  that: 

at  the  lower  wall: 
urij_!/2,j  =<Jr,j+ 1/2,1 

Pj- l/2,i  =Pj+l/2,i  (4.4.1-la) 

Pj-i/2,i  =Pj*  1/2,1 

vV  1/2,1  1/2.1 

Y\- !  /2,i  =-yj  +  ]/2>] 

t|j-l/2,i  1/2,1 

wherefore 

^>i/2,i  =sign'6~j+  j/2,j 

^H/2,1  =sign‘5-j-1/2,i  (4.4.i-2a) 

at  the  upper  wall: 

u/ij+,/2ij  =u^_  1/2,1 

Pj+l/2,i  =Pj-!/2,i  (4.4.1-lb' 

Prl/2,i  =Pj-i/2,i 

vnj+1/2,i="^nj-l/2,i 


wherefore 


=sign"5'j-i/2ti 

=sign'6-j+i/2ij 


(4.4. 1 ->D) 


ri 


"5>  1/2,1 


where 


siqn=  -1  for  the  mass,  x-momentum  and  energy  fluxes. 
♦  1  for  the  y-momentum  fiux. 


in  the  case  of  the  H-grld,  on  the  other  hand,  the  same  simple 
implementation  of  a  mirror  image  at  the  boundary  is  clearly  not  correct. 
This  can  for  example  be  seen  in  Fig.  4. 4.1-1,  where  a  discontinuity  in  the 
grid  is  introduced  at  the  boundary  for  the  mirror  image,  it  is  seen  by 
considering  the  elements  created  by  the  mid  grid  points  throughout  the 
computational  domain  that,  at  the  boundary,  the  flow  which  should  be 
calculated  in  point  (j,i)  is  instead  calculated  for  a  point,  in  tne  present 
example,  at  (j,i+5). 

For  the  H-grid,  the  boundary  conditions  are  insteac  implemented  in 
different  ways  for  the  option  that  the  grid  points  are  situated  directly  at 
the  computational  boundaries.  The  F’  flux  can  be  calculated  as  in  the 
interior  points,  as  no  information  from  outside  the  computational  domain  is 
necessary  for  this  flux  at  the  blade  walls  (obviously,  at  the  inlet  and 
outlet  F"  would  need  information  from  outside  the  computational  domain, 
whereas  G‘  would  not  need  any  extra  information).  For  determining  tne  G‘ 
flux,  different  possibilities  exist.  One  (Method  1)  would  be  to  extrapolate, 
linearly  or  parabolicly,  the  TW“  fluxes,  from  the  splitted  fluxes  inside  the 
computational  domain,  towards  the  first  mid  point  outside  the 
computational  domain.  This  would  be  consistent  with  the  way  the  fluxes 
are  determined  in  the  rest  of  the  computational  domain,  and  would 
introduce  a  numerical  dissipation  at  the  walls  similar  to  the  rest  of  the 
flow  field.  Another,  sligthly  different  approach  (Method  2),  consists 
instead  of  extrapolating  the  splitted  fluxes  <W~  from  the  mid  points,  we 
instead  compute  the  total  fluxes  in  the  gridpoints  close  to  the  walls  and 
extrapolate  G\  in  the  first  mid  grid  point  outside  the  computational 
domain,  parabolicly  from  these.  This  gives: 


For  the  upper  wall: 


>V,Vy. 


GVl/2.r(36’j-2-r  10-G*j_i  ,i4  156*Jij)/8 


(44.1 -3a) 


For  the  lower  wall: 


G  i-i/2, i=^6  j»2,rlOG  ]+i,i+^6  j j)/8 


(4.4! -3b) 


Furthermore,  in  this  approach,  we  are  also  missing  the  information  of 
(5*j+  i/2(j  at  the  lower  wall  calculation  and  the  T5~j-i/2.i  at  the  upper  wai! 
calculation.  This  is  seen  from  the  way  the  '5±j±i/2ii  fluxes  are  calculated 
with  the  second  order  spatial  accuracy  (see  section  3.2).  Therefore,  for 
calculation  of  the  G’n  at  the  walls  and  the  point  next  to  the  wall,  the  total 
flux  in  point  (i*  1/2,0,  for  the  lower  wall,  and  point  (J- 1/2,0,  for  the  upper 
wall,  is  computed  as: 


Lower  wall:  G'j+|/2ij=(G‘jj+G*j.*.|i|)/2 
Upper  wall:  G'j_|/2ii=(G*j  i+G‘j_]  \)/2 


(44.1-4) 


respectively1. 

Obviously,  the  so  computed  flow  variables  at  the  walls  do  not  satisfy  the 
boundary  condition  of  no  mass  flow  through  the  wall.  Therefore,  the 
Riemann  invariant  normal  to  the  wall  (R+  at  the  upper  wall,  R~  at  the 
lower  wall,  see  for  example  Appendix  A4)  is  calculated  from  the  computed 
flow  variables.  The  introduction  of  the  boundary  condition  (zero  normal 
velocity  at  the  wall)  is  thereafter  used  together  with  the  Riemann 
invariant  to  determine  the  relationship  between  normal  velocity  3nd 
pressure  at  the  wall. 

These  methods  of  extrapolation  for  the  boundary  is  known  from  the 


author's  previous  work  (with  a  MacCormack  predictor-corrector 


without  explicit  numerical  dissipation)  111  to  give  good  results  as  regards 
to  the  pressure  and  density.  Howeve. ,  some  error  can  be  introduced  (with 
tne  magnitude  depending  on  the  angle  between  the  £=const  ana  q=consr 
lines)  in  the  stagnation  pressure.  Therefore,  the  velocity  tangential  to  me 


1  Such  3n  averaging  must  be  carefully  checked  as  it  eventually  may 
introduce  errors  and/or  instabilities  in  some  cases. 


blade  surface  can  be  updated,  in  steady  state  flow,  by  considering  that  the 
wall  corresponds  to  a  streamline,  wherefore  tne  entnaipny  is  Known  aiong 
it.  This  gives  the  following  relationship,  as  v=0  (velocity  component 
normal  to  the  wall): 

u2=2{hc_00-a2/[^f- 1  ]}  (4.4. 1-5) 

with  u=  velocity  component  tangential  to  the  wall. 

However,  previous  work  of  extrapolation  at  the  boundaries  has  not  included 
splitting  of  the  fluxes,  and  thus  hot  any  special  treatment  for  the  points 
(j+l/2,i)  for  the  lower  wall  calculation  and  the  point  (j- 1/2,0  for  tne 
upper  wall  calculation.  Simple  averaging  of  the  total  fluxes  as  in  eq. 
(4.4. 1-4)  may  eventually  lead  to  instability. 

The  experience  with  the  present  flux  vector  splitting  method  has  shown 
that  the  Method  2  of  implementing  the  boundary  conditions  gave  good 
results  for  the  first  two  dimensional  test  problem  (section  4.6.1),  both  for 
the  O-grid  and  tne  H-grid  (at  least  as  long  as  the  boundary  conditions  at 
the  inlet  and  outlet  are  overspecified),  but  that  most  of  the  cases  for  the 
second  test  problem  (section  4.6.2)  did  not  converge. 

4.4. 1. b.  Grid  Points  Situated  1/2  Point  from  the  Walls. 

To  eliminate  these  convergence  problems,  it  was  instead  tried,  with  better 
resuits,  to  implement  the  flux  vector  splitting  methodology  in  a  finite 
volume  direction,  althougn  the  actual  solution  aigontnm  is  stiii  a  finite 
difference  approach. 

The  grid  was  then  instead  constructed  such  that  the  cell  interface  lies  on 
the  boundary  of  the  computational  domain  (section  4.2.2). 

The  boundary  conditions  are  then  incorporated  as  already  explained  in 
section  4.2.2. 

4.4.2.  Boundary  Conditions  at  the  Inlet. 

This  is  performed  in  exactly  the  same  way  as  in  one  dimensions  (section 

3.3.2. P)  with  the  additional  constraint  of  the  inlet  flow  direction. 

4.4.J.  Boundary  Conditions  at  the  Outlet. 


As  in  one  dimensional  flow  (section  3. 3. 2. a)  with  the  additional  condition 
tnat  v  is  accepted  from  the  calculation,  or  tas  a  second,  non-physicai, 
option)  that  v=0. 
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Fig.  4.4. 1- 1 .  illustration  of  inaccuracy  of  mirror  linage 

treatment  of  tne  oounoary  for  an  H-gna. 
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4.5.  Numerical  Algorithm  In  Two  Dimensions. 

The  solution  algorithm  is  the  same  as  in  one-dimensional  flow,  with  the 
second  dinTension  added. 


First  step: 


WjjK  -At{A(Fj  j+)k+A(Fj  j-)k}/Ax 
-At{A(Gjij+)^A(Gjif)k)/Ay 


Second  step: 


wj.r1 


-At[A(F*jj+)k+A(F,‘j|j“)k]/Ax 
-At[A(G*  j  j*  )k*  A(G**  j  j-)kJ/Ay} 


with  the  spatial  differencing  (in  both  dimensions,  one  at  a  time)  as  in 
section  3.4. 


I_ 


.  ....... 
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4.6.  i .  First  Test  Problem. 


The  first -two  dimensional  test  problem  is  identical  to  tne  tniro  one 
dimensional  problem,  i.e.  a  duct  with  radial  subsonic  flow  (Fig.  4.6.1- 1) 
with  an  outlet  Mach  number  of,  at  the  centerline,  M=0.8864.  For  the 
purpose  of  a  detailed  investigation  of  the  influence  of  the  number  of  grid 
elements  on  the  numerical  solution,  the  opening  of  the  duct  was  chooser,  to 
be  fairly  large  (wedge  angle=90°).  Furthermore,  two  sets  of  grids  were 
investigated.  Tne  first  is  a  O-grid  (Fig.  4.6.1 -la)  and  the  second  an  H-gno 
(Fig.  4.6.1 -lb).  As  tne  flow  in  the  case  under  study  is  radial,  it  is 
obviously  expected  that  the  first  grid  will  give  the  most  accurate  resu't 
for  a  predefined  number  of  grid  points. 

The  results  from  the  two  dimensional  O-grid  is  illustrated  in  Fig.  4.6. 1 -2. 
Here,  the  lines  of  constant  Mach  number  are  compared  to  the  exact  results. 
(For  the  first  test  problem  all  computations  shown  are  performed  with 
overspecified  boundary  conditions  at  the  inlet  and  outlet.)  The  comparison 
is  done  in  the  physical  plane  for  a  general  view  (Fig.  4.6.1  -2a)  ano  in  tne 
computational  plane  for  more  details  (Fig.  4.6.1 -2b).  The  symmetry  of  tne 
solution  is  exact.  However,  as  already  concluded  for  the  quasi  ore 
dimensional  computations,  21  grid  points  in  the  streamwise  direction  ’$ 
not  enough  to  accurately  reproduce  the  exact  solution  (due  to  the  numerical 
dissipation  inherent  in  the  flux  vector  splitting).  With  21  streamwise 
points  (and  9  circumferential)  the  numerical  solution  is  shifted 
approximately  I  element  at  Isomachline  M=0.5  (a  similar  shift  is  present 
also  in  the  density  and  pressure).  Furthermore,  as  expected,  3  grid  points 
(=2  elements)  in  the  circumferential  direction  gives  a  very  inaccurate 
solution  (Fig.  4.6.1 -2b),  whereas  9  and  13  circumferential  points  give 
similar  results. 

For  the  H-grid,  results  are  presented  in  Fig.  4.6. 1-3  for  the  case  with  tke 
boundary  conditions  at  the  walls  implemented  according  to  Method  2 
section  4.4. i. a  (extrapolation  of  total  fluxes  to  a  midpoint  outside  the 
computational  domain).  For  the  computational  grid  consisting  of  60x i 2 
elements  there  is  virtually  no  difference  between  the  exact  and  comouteo 
solutions,  if  the  tangential  velocity  component  at  the  wall  is  corrected  as 
part  of  the  boundary  condition  implementation  according  to  eq.  (4.4.1 -5). 
This  is  shown  for  the  density  (Fig.  4.6. 1 -3a),  pressure  (Fig.  4.6.1 -3b)  and 
Mach  number  (Fig.  4.6.1 -3c).  if  the  correction  of  u  is  not  performed,  a 
slight  difference  between  tne  computed  and  exact  results  is  found  :n  the 
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upper  left  corner  of  the  nozzle  for  the  density  and  pressure  (Fig.  4.6.1- 
3a,D).  The  isomachlines  show  however  in  this  case  a  larger  deviation  from 
the  exact  values  (Fig.  4.6.1 -3c).  This  is  due  to  the  fact  that  the  stagnation 
pressure  is 'overestimated  at  the  wall.  This  is  a  local  effect,  and  the 
turning  of  the  Isomachlines  is  limited  to  the  gridpoints  at  the  wall.  This 
is  seen  from  Fig.  4.6.1 -3c,  where  the  values  agree  in  all  points  inside  the 
flow  field. 

A  stability  factor  of  STAB=0.4  (see  eq.  6.0-4,  6.0-6)  was  used  for  all 
these  calculations. 
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4.6.2.  Second  Test  Problem. 

The  second  two  dimensional  test  problem  is  the  one,  by  now  standard  (/32- 
35/,/7/),  consisting  of  a  straight  channel,  with  a  l  OX  double  circular  arc 
bump  on  one  of  the  walls  (Fig.  4.6.2- 1).  Computed  lnviscld  results  with 
this  geometry  have  been  presented  by  several  authors,  for  different  inlet 
flow  conditions,  since  its  introduction  by  6.  Morettl  in  1981  /3 2/.  No 
experimental  data  can  be  compared  to  results  from  the  invlscid  methods  as 
the  viscous  flow  separates  at  about  70X  of  the  bump  length  for  an  inlet 
Mach  number  of  M.^0.5  /36/. 

Calculations  with  the  present  method  where  performed  with  an  H-grid  at 
inlet  Mach  numbers  of  M-^.S  (completely  subsonic  flow)  and  M-^0.675 
(transonic  flow). 

The  first  case  (M.^0.5)  was  computed  using  only  6x20  elements  (7x21 
grid  points)  and  with 

grid  lines  */the  walls  (section  4.2.1) 

boundary  treatment  of  r)-f luxes  according  to  method  2  in  section 
4.4.1. a. 

This  method  gave  results  (for  overspecified  boundary  conditions  at  inlet 
and  outlet)  in  good  agreement  with  the  ones  presented  in  /7,33/,  but 
Instabilities  was  found  for  other  thicknesses  of  the  bump,  when  the 
overspecification  of  the  boundary  conditions  at  the  inlet  and  outlet  was 
taken  away. 

The  second  case  (M.^sO.675)  was  then  instead  computed  with  the  grid 
generated  so  that  the  wall  Is  located  t/2  cell  away  from  the  grid 
points  { section  4.2.2),  with  corresponding  boundary  conditions  at  the 
walls  (section  4.4.i.b). 

The  pressure  on  the  wall  with  the  bump  is  presented  in  Fig.  4.6.2-2a.  For 
this  computation  the  pressure  at  the  walls  was  calculated  accoiding  to  the 
zeroth  order  extrapolation  in  section  4.4. l.b,  t.e.  Pj*i/2,i=Pj.i  with  jti/2 
located  on  the  upper  and  lower  walls,  respectively.  The  first  calculation 
(Run  201)  was  performed  with  overspecified  boundary  conditions  at  inlet 
and  outlet,  and  the  second  (Run  202)  with  radiative  boundary  condition  at 
inlet  and  constant  pressure  at  outlet  (implemented  with  "Method  1"  in 
section  3.3.1).  Both  computations  are  performed  without  special  treatment 
at  the  first  point  downstream  of  the  shock,  l.e.  second  order  spatial 
accuracy  (flux  limiter  0=  1 ,  see  section  3.2)  everywhere,  and  with  the 


entropy  in  the  two  points  closest  to  the  wall  put  equal  to  the  entropy  in 
the  third  point  away  from  the  wall  in  order  to  reduce  numerical  entropy 
errors  at  the  wall  (see  for  example  /11;37/).  There  is  virtually  no 
difference  between  the  solutions,  apart  from  close  to  the  inlet  and  outlet 
boundaries. 

The  highest  Mach  number  is  approximately  Mm^l.38  and  is  reached  at  70% 
of  the  bump  length.  This  is  In  close  agreement  with  results  from  /33/ 
1  -38,  72%)  and  /34/  (Mmax*1-32,  72%),  whereas  there  is  some 
difference  towards  the  results  presented  in  /35/  (Mmax^l.5,  70%). 

In  Fig.  4.6.2-2b  the  lines  of  constant  isentropic  Mach  numbers  are 
presented  In  the  computational  domain  for  the  same  calculations.  The 
results  are  again  in  agreement  with  ref.  /33,34/. 

It  should  however  be  noted  that  the  results  introduce  an  Increase  of 
stagnation  entalphy  at  the  wall  with  the  bump  in  the  case  of  the  radiative 
boundary  condition  at  the  inlet  (compare  /24/  where  this  also  occur).  This 
increase  corresponds  to  Ahc*2%,  with  an  increase  also  at  the  inlet.  If  hc  is 
instead  specified  at  the  Inlet  (for  example  as  In  the  “capacity"  boundary 
condition)  the  level  of  error  in  enthaalpy  stays  approximative^  the  same, 
but  there  is  no  increase  of  hc  at  the  inlet. 

if  the  stagnation  enthalpy  at  inlet  and  outlet  are  not  the  same  due  to 
numerical  errors,  the  Mach  number  as  calculated  will  also  have  numerical 
errors  (compare  for  example  the  previous  discussion  on  this  matter  in 
section  3.3.1.b),  even  if  the  pressure  is  correct.  More  points,  and  clustered, 
close  to  the  wall  would  improve  this  inaccuracy,  as  would  a  better 
extrapolation  of  the  pressure  at  the  wall  (section  4.2.2). 

Furthermore,  a  grid  which  Is  orthogonal  to  the  walls  would  increase  the 
accuracy. 
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Lines  of  constant  Mach  number 
(determined  from  the  pressure) 
in  the  comDutaJsenal  plane. 
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Fig.  4.6.2-2.b.  Second  two  dimensional  test  problem:  Lines  of 
constant  Mach  number  (Grid  15x6 1 ). 
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Part  II:  Testing  of  the  Methodology  for 
Unsteady  Flows  in  Non-Moving  Grids. 


.ny* 


5.  One  Dimensional  Unsteady  Test  Problems  and  Computational 
Results. 

Before  going  into  the  application  of  the  unsteady  flow  in  a  moving  cascade 
the  algorithm  should  first  be  tested  on  unsteady  flows  In  a  fixed  reference 
system.  For  these  tests,  two  cases  were  chosen.  The  first  (section  5.1)  is 
the  shock  tube  problem  in  order  to  determine  if  the  shock  speed  as 
captured  by  the  method  corresponds  to  the  exact  solution,  and  how 
accurately  the  moving  shock  and  expansion  fans  can  be  determined.  The 
second  (section  5.2),  concerns  the  unsteady  shock  movement  in  a  quasi  one 
dimensional  nozzle  with  oscillating  back  pressure. 

It  should  be  pointed  out  that  all  results  presented  in  this  section  are  based 
on  the  original  flux  vector  splitting  by  Steger/Warming,  with  the  spatial 
differencing  performed  according  to  the  "MUSCL"  approach  (see  sections  3 
and  4  for  details) 


5.1.  Shock  Tube. 


The  shock  tube  chosen  corresponds  to  the  one  employed  in  /l 3/.  The  ratio 
of  the  high  pressure  (left  region  in  Fig.  5.1-1)  and  low  pressure  (right 
region)  is  10.  The  diaphragm  is  situated  at  x=0.95  and  bursted  at  time 
t=0.  The  solution  of  density  at  a  time  t=!  is  given  in  Fig.  5.1-1.  The 
following  information  can  be  drawn: 

The  positions  of  the  shock,  contact  discontinuity  and  expansion  fan 
are  accurately  captured,  which  means  that  the  speeds  of  these  waves  are 
correctly  calculated. 

The  expansion  fan  starts  with  an  overshoot  of  approximately  4%  and 
is  smothed  out  over  approximately  15  points  instead  of  12  points  for  the 
exact  solution. 

The  contact  discontinuity  is  completely  smeared  out  and  extends 
over  6  points. 

The  shock  is  accurately  captured  with  a  sllgth  undershoot  before 
and  after,  and  It  extends  over  two  elements  only  (one  point  in  the  shock 
region). 

The  density  in  all  four  regions,  separated  by  the  different  surfaces, 
is  accurately  predicted. 

It  can  thus  be  concluded  that  the  present  methodology  can  predict  the 
correct  shock  speed  and  density  ratio  in  a  shock  tube  problem. 


\ 


5.2.  Shock  Oscillation  Induced  By  Time  Varying  Back  Pressure. 

5.2. 1 .  validation  of  tne  netnoq. 

Adamson  et  al  have  developed  analytical  small  perturbation  theories  to 
the  problem  of  oscillated  back  pressure  In  transonic  one  3r.d  two 
dimensional  nozzles  /39-45/.  Their  results  indicate  that  the  unsteady 
shock  oscillation  amplitude,  as  well  as  the  time-averaged  shock  position, 
depends  not  only  on  the  steady-state  (or  average  back  pressure)  but  also 
on  the  amplitude  of  the  oscillation  and,  largely,  on  the  oscillation 
frequency  (and,  of  course,  the  nozzle  geometry).  As  these  solutions  are  of 
an  analytical  nature,  comparisons  with  the  present,  fully  numerical, 
methodology  is  of  interest. 

The  Adamson  et  al  nozzles  are  prescribed  by  the  following  function  of  the 
area  /44/: 


A(x)=W2.f(x) 


where  f(x)  is  a  function  describing  the  wall  shape  and  e  is  a  constant 
defining  the  small  perturbation.  In  ref.  /44/  the  function  f(x)  is  defined 
as,  with  the  entrance  of  the  subsonic  part  of  the  nozzle  situated  at  x=- 
0.5,  the  throat  at  x=+0.5  and  the  exit  at  x=1 .5: 


3 

-0.5 

SX 

S- 1/6 

{27[3(x-0.5)+2P-48[3(x-0.5)+2P}/ 1 3*3 

-1/6 

sX 

<1/6 

18[3(x-0.5)]2/l3 

1/6 

iX 

s5/6 

(27(3(x-0.5)-2P+48[3(x-0.5)-2P)/  1 3*3 

5/6 

sX 

s7/6 

3 

7/6 

SX 

<1.5 

Furthermore,  the  unsteady  flow  is  introduced  by  varying  the  back  pressure 
according  to 


p=pe+peiampSin(w[t-to]) 


with 


w=circular  frequency  of  oscillation=2n/tcycie 


The  first  results  presented  compare  the  numerically  computed  shock 
position  with  the  analytical  location  from  ref.  /44/,  with  the  following 
(dimensionless)  parameters: 

£  =0.1 

pe  =0.626474 

Pe,a!Tip=Q-0 14032  (i.e.  Pe,amp/Pe=2-2%) 
tcycle  =  ,S  2  =>  k  =Coolr]/(2qrl 

=[27T]/[2tcycie)‘V/(trQrl 

={2TT]/(2t(;yc1e] 

=0.173 

(if  the  following  reference  values  are  employed: 

Tr  =300K 

lr  reference  length 

qr  reference  velocity 

tr  =lr/Qr 

we  find  the  physical  values  to  be: 

tcycle  .physical  =tcycle’tr  *0.006207  S 

f  =  1  /  ^cycle  *161  Hz 

Based  on  a  reference  length  of  !r=lm,  the  frequency  would  instead  be 
f=16. 1  Hz.) 

In  Fig.  5.2- la  the  numerical  (121  grid  points)  and  exact  steady  state 
pressure  distributions  are  given.  As  for  the  previous  steady  state 
solutions  (section  3.5)  the  shock  is  captured  sharply  and  accurately,  if  tne 
position  of  the  sonic  line  within  the  structure  of  the  numerical  shock  is 
defined  as  the  shock  location  (as  pointed  out  in  /43/  to  be  a  reasonable 
assumption),  It  is  concluded  that  the  exact  and  numerical  shock  positions 
agree  well. 

As  the  backpressure  then  starts  to  oscillate  (at  a  dimensionless  time  of 
t*4i . I ,  Fig.  5.2- 1b),  the  shock  will  not  react  until  the  disturbances  wave 
has  travelled  the  distance  from  the  outlet  (x=i.5)  to  the  shock  (x=  l .00). 
Tne  snock  tnen  starts  to  react,  in  the  beginning  slowly,  to  tne  change  in 


=0. 1  m 

=/RT  *293  m/s 
*3.41  -JO-4  s 


pressure  at  an  approximate  time  of  t*  44  (Fig,  5.2- 1  c).  As  the  pressure 
increases  in  the  beginning  of  the  cycle  (Fig.5.2-1b)  the  shock  moves,  as 
expected,  upstream,  it  is  seen  that  the  shock  position  as  determined  by 
the  numerical  method  agrees  well  with  the  analytical  small  perturbation 
from  Adamson  et  al  /44/,  although  the  numerical  solution  displaces  the 
average  shock  location  sllgthly  downstream  (Axs*0.0l2,  corresponding  tc 
*0.7  mesh  elements  with  the  121  grid  points  used). 

Both  methods  illustrate  the  important  factor  that  the  average  shock 
position  of  an  unsteady  shock  does  not  have  to  be  the  steady  state 
position.  In  the  example  presented  in  Fig.  5.2- 1  c  the  steady  state  and  the 
averaged  unsteady  (based  on  the  first  cycle  only)  shock  positions  are  as 
follows: 


Analytical1 

Numerical 


Steady  state 
shock  position 
1.000 
0.998 


Unsteady  average 
shock  position 
0960 
0.972 


The  corresponding  Mach  number  distributions  for  different  time  levels 
during  the  first  cycle  are  given  In  Fig.  5.2-1d.  This  figure  illustrates  the 
shock  displacement  in  the  channel.  It  is  also  seen  that  the  flow  upstream 
of  the  shock  is  not,  as  expected,  Influenced  by  the  back  pressure 
oscillations.  Also  in  this  figure  it  is  clearly  seen  that  the  average  shock 
location  is  not  identical  to  the  steady  state  position. 

The  reason  for  this  is  probably  the  variation  in  the  flow  velocity  relative 
to  the  shock,  as  well  as  the  phase  lag  between  the  post  shock  pressure 
change  and  the  shock  position,  combined  with  the  pressure  fluctuation 
downstream  of  the  shock.  As  the  shock  moves  upstream  in  the  sample 
case  shown,  the  flow  velocity  relative  to  the  shock  Increases  (as  for  this 
sample  case  the  shock  velocity  becomes  larger  than  the  change  in  the 
velocity  upstream  of  the  shock  due  to  the  area  change  (see  ref.  /44/  as 
well  as  the  later  discussion  on  shock  velocity)).  This  increases  the 
pressure  jump  over  the  shock,  at  the  same  time  as  the  disturbance  wave 
tram  from  the  exit  introduces  a  compression  In  the  post  shock  flow  field, 


1  Please  not  that  it  Is  not  entirely  clear  from  ref.  /44/  if  the 
analytical  solution  was  performed  over  several  cycles  or  if  it  was 
stopped  when  the  pressure  returned  to  the  steady  state  value. 


slightly  (as  the  change  in  nozzle  area  is  small)  counteracted  by  an 
expansion  due  to  the  increase  in  the  nozzle  section  up  to  x=7/6  (Fig.  5.2- 
id,e>. 

After  the  shock  has  reached  its  minimum  position  (x$imin*0.930),  it  starts 
to  move  downstream  again.  At  a  certain  time  it  passes  again  over  its 
unsteady-  average  location  (xs, unsteady  average*0-972,  see  Fig.  5.2- ic).  In 
this  position  the  shock  has  reached  its  maximum  speed  going  downstream, 
wnerefore  tne  flow  velocity  relative  to  the  shock  has  its  minimum  value 
Simultaneously,  tne  disturbance  wave  train  from  the  nozzle  exit 
introduces  an  expansion  downstream  of  the  shock.  The  pressure  jump  at 
the  shock  is  thus  smaller  than  the  jump  when  the  shock  moved  upstream 
(Fig.  5.2-1  e). 

With  these  two  phenomena  acting  together  (in  different  ways,  depending 
on  the  nozzle  geometry,  pressure  amplitude  and  frequency)  it  is  not  to  be 
expected  that  the  average  of  the  unsteady  shock  position  should  coincide 
with  tne  steady  state  shock  location. 

if  the  numerical  calculation  is  continued  over  several  periods  it  is  seen 
that  a  periodic  shock  movement  is  found  after  approximative^  two  cycles 
(Fig.  5.2- if). 

For  a  lower  frequency  of  the  Imposed  exit  pressure  oscillation,  but  with 
the  same  amplitude  as  previous  (f=32Hz,  k=0.035  based  on,  as  before, 
lr=0.1m),  the  shock  movement  increases  drastically  (compare  Figs.  5.2- If 
and  5.2-23).  This  can  be  explained  by  the  fact  that,  for  high  reduced 
frequencies,  the  shock  has  not  got  time  to  fully  react  to  the  disturbances 
(for  example  a  compression  wave),  before  there  is  another  disturbance 
(for  example  an  expansion  wave)  counteracting  the  effect  of  the  first  one. 
As  in  the  previous  sample  case,  It  is  determined  that  the  unsteady  flow  is 
practically  periodic  after  two  cycles.  Furthermore,  the  unsteady  average 
value  of  the  shock  position  has  again  changed,  and  the  following  values 
are  obtained: 


Average  shock  position 

(numerical, 

based  on  first  cycle) 

0.998 


1.013 
0  972 


4x$= 

x$, unsteady'  xs,  steady 


Steady  state 
k=0.035 
k=0  I  73 


0.015 

-0.026 


The  increase  in  shock  movement  can  also  be  seen  by  studying  the  Mach 
number  distribution  along  the  nozzle  at  the  steady  state,  maximum  and 
minimum  positions  for  the  first  cycle  for  the  two  sample  cases  (Figs. 
5.2- Id,  5.2-2b).  In  the  low  reduced  frequency  case  the  maximum  shock 
position  extends  well  into  the  part  of  the  nozzle  where  the  area  is 
constant  (x*7/6).  The  minimum  and  maximum  Mach  numbers  are 
approximatively  1. 1 1  and  1.20  resp.,  compared  to  the  steady  state  value 
of  1.18,  and  the  shock  moves  over  approximatively  36%  of  the  nozzle 
distance  downstream  of  the  throat,  compared  to  8%  for  the  higher  reduced 
frequency. 

in  both  cases  the  imposed  pressure  oscillation  is  2.2%. 

Also  for  the  lower  reduced  frequency  the  unsteady  shock  position,  as 
calculated  with  the  flux  vector  splitting  method,  agrees  well  with 
analytical  results  from  Adamson  and  Liou  /44/  (see  Fig.  5.2-2c). 

As  a  final  validation,  and  before  applying  the  flux  vector  splitting  method 
to  a  few  different  one  dimensional  sample  cases,  the  shock  velocity  is 
compared  with  the  analytical  results  in  /44/.  it  should  here  be  noted  that 
the  numerically  determined  shock  velocity  is  the  time  derivative, 
calculated  as  a  centered  difference,  of  the  shock  location,  i.e. 

us  =  Axs/At 

Furthermore,  xs  is  an  interpolated  value,  given  by  the  position  of  the  sonic 
line  within  the  numerical  shock.  It  is  thus  not  expected  that  the 
numerical  shock  velocity  should  be  uniform,  but  rather  that  us 
fundamental  frequency  should  coincide  with  the  analytical  shock  velocity 
from  /44/, 

This  comparison  is  made  in  Figs.  5.2-3a,b  for  k=0. 173  and  k=C.C35, 
respectively.  In  both  figures  the  curve  from  Adamson  and  Liou  /44 /  agrees 
well  with  the  average  of  the  numerical  results.  However,  the  bandwidth  of 
the  numerical  solution  is  considerable,  which  can  be  explained,  in  a 
physical  manner  considering  the  numerical  scheme,  by  looking  at  Fig  5  3- 
2a. 

It  must  be  recognized  that,  although  the  numerical  shock  is  fairly  sharp, 
it  still  consislts  of  a  few  intermediate  grid  points.  The  disturbances 
moving  upstream  from  the  exit  will,  as  they  approach  the  two  or  three 


grid  points  immediately  downstream  of  the  shock  position  xs  (which  is, 
we  repeat,  determined  as  the  position  of  the  some  line  within  the 
numerical  shock  structure),  find  a  higher  flow  velocity  (see  steady  state 
solution,  Fig.  5.2- la).  The  speed  of  propagation  of  the  disturbance  signals 
(=u-a)  will  thus  slow  down.  As  the  last  point  upstream  of  the  shock 
starts  to  be  influenced,  the  shock  velocity  rapidly  increases  (at  t~44.  l), 
until  t*45.7.  At  t~46.4  the  Mach  number  in  grid  point  1=90  passes  from 
supersonic  to  subsonic.  The  pressure  in  grid  point  1=89  starts  then  to 
build  up  (with  an  initial  low  shock  velocity)  until,  at  t*47.9,  the  Nacn 
number  in  this  grid  point  becomes  subsonic.  At  this  time,  the  flow  in  grid 
point  1=88  starts  to  be  influenced,  etc.  A  close  look  on  the  shock  position 
curve  in  Fig.  5.2- 1c  reveals  also  this  jump  every  time  the  Mach  number 
passes  the  sonic  value  in  a  grid  point. 

The  shock  velocity  traces  at  the  lower  reduced  frequency  (k=0.035,  Fig. 
5.2-3b)  does  not  show  an  identical  behaviour.  However,  this  is  just  due  to 
the  fact  that  only  each  i  00th  time  step  was  used  to  compute  the  shock 
velocity  (and  also  the  shock  position  in  Fig  5.2-2a) 

Another  important  aspect  of  the  oscillating  shock  is  when  it  reacts  as 
regards  to  the  imposed  pressure  fluctuation,  it  has  been  shown 
analytically  (for  small  perturbation  pressure)  by  Adamson  and  Llou  /4 4/ 
that  the  two  sample  cases  presented  previously  corresponds  to  two 
d’fferent  classes  of  solutions: 

The  first  corresponds  to  a  case  when  the  disturbance  wave  from  the 
exit  reaches  the  shock  in  a  time  comparable  to  the  time  the  imposed 
pressure  performs  a  cycle 

The  second  corresponds  to  the  case  when  the  period  of  the  imposed 
pressure  is  much  longer  than  the  time  it  takes  for  the  disturbance  wave 
to  reach  the  shock? 

Each  of  thesr  two  classes  has  a  different  analytical  solution  (see  for  ex. 
ref.  /44/,  p.25),  with  the  implication  that,  for  the  first  solution  (i.e. 
tcycie/tsh*^1^.  the  shock  velocity  is  essentially  in  phase  (apart  from 
the  time  'ag  ts^)  with  the  imposed  pressure  fluctuation,  whereas  for  the 

*  irTthe  first  sample  case  the  time  for  the  disturbance  wave  to  reach 
the  snock  is  approx im at ively  ^*=43  4-41  1*2.3  (see  Fig  5.2-1  c)  ana  the 
penond  tCyCie=18.2,  i.e.  tCyCie/tsh*8.  In  the  second  sample  case,  t5fl  is 
still  the  same,  whereas  tcycie*91.  Thus,  tcycie/tsri*40. 


second  solution  (i.e.  tcyC]e/tSh»  1 )  the  phase  lag  depends  upon  other  terms, 
such  as  the  nozzle  geometry  and  average  exit  pressure. 

This  implies  that,  for  the  case  presented  in  Fig.  5.2-1,  the  shock. 
position3  should  be  displaced  as  regards  to  the  pressure  fluctuation  at 
the  outlet  approximative^  At*tcyC)e/4+tsh»7.  As  the  maximum  pressure 
(Pe,max>  appears  at  tpe.max*^.?  (Fig.  5.2- lb),  the  minimum  shock,  position 
snould  be  reached  at  tmin  sh0C|Ca*tpe>max*At«52.7.  This  is  the  case 
(obviously  for  the  analytical  solution,  but  also  for  the  numerical),  as  can 
be  seen  from  Fig.  5.2- 1c. 

However,  as  explained  in  /44/,  for  the  second  sample  case  (Fig.  5.2-2) 
another  phase  lag  should,  and  does  Indeed,  exist.  The  analytical  solution 
by  Adamson  and  Liou  (ref.  /44/,  p.29)  shows  that  the  shock  velocity  for 
iow  reduced  frequencies  (tcycie/tSh»l)  does  not  respond  instantenuous  to 
a  pressure  disturbance,  although  the  signals  from  the  imposed  pressure 
fluctuation  at  the  exit  reaches  the  shock  “  instantenuous  ly“  (as  tSh  now  is 
very  small  compared  to  our  time  scale  of  interest,  i.e.  the  period  of 
oscillation).  This  is  due  to  some  terms  in  the  small  perturbation  analysis, 
which  are  dependent  on  the  nozzle  form,  pressure  amplitude  and 
frequency,  etc.  /44/.  it  is  important  to  note  that  the  numerical  solution 
accurately  predicts  this  phase  lag,  together  with  the  amplitude  of  the 
shock  oscillation  (Fig.  5.2-2.C). 
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Note:  p*sin(t)=>us*sln(t)=>xs*cos(t)  for  this  case. 
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as  seen  from  the  two  comparisons  with  ref.  / 44/  the  shock  may  move 
considerably,  even  if  the  imposed  pressure  fluctuation  at  the  exit  is 
relatively  small,  in  such  a  case,  for  a  given  nozzle  geometry,  the  shock 
movement  3nd,  in  some  cases  as  will  be  seen  later,  the  flow  throughout 
the  entire  nozzle  depends  on  the  reduced  frequency  of  the  prescribed 
disturbance.  In  the  following  we  will,  for  the  same  nozzle  geometry  as 
previously,  investigate  the  flow  field  throughout  the  channel  for  a  higher 
back  pressure  fluctuation  in  three  sample  cases: 


Sample  peH 
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Sample  Case  J: 

As  previous,  the  imposed  pressure  fluctuation  is  sinusoidal,  with  an 
initial  increase  at  t=t0  5.2-4. 1).  The  pressure  fluctuation  is  such 
tnat  if  a  steady  state  pressure  was  kept  at  the  maximum  exit  pressure, 
the  fiow  would  have  been  subsonic  in  the  whole  nozzle. 

The  change  in  pressure  distribution  in  the  nozzle  throughout  the 
oscillation  is. shown  in  the  sequence  of  Figs.  S.2-4.2  to  5.2-4.32,  with  12 
plots  per  period  of  oscillation.  Initially  the  imposed  fluctuation  at  the 
exit  travels  upstream  (Fig.  5. 2-4.3).  The  pressure  downstream  of  the 
shock  is  thus  modified,  essentially  according  to  the  imposed  pressure 
oscillation  (since  the  nozzle  has  a  constant  cross  section  for  x>7/6).  At 
the  approximate  time  t*n/2  (where  the  time  has  been  referenced  to  the 
start  of  the  pressure  oscillation  for  easier  identification),  the  wave  tram 
has  reached  the  shock  (Fig.  5. 2-4.5)  and  has  influenced  all  the  post  shock 
computational  points  (as  well  as  the  first  supersonic  point,  which  is 
situated  within  the  numerical  shock  structure).  At  the  next  time  level 
presented  (t=4n/6)  the  post  shock  pressure  has  increased  as  the.  wave 


train  has  propagated  upstream  (Fig.  5.2-46).  The  shock  has  started  to 
move  slightly  upstream.  At  t«5n/6  the  pressure  just  downstream  of  the 
shock  has  reached  its  maximum  level  (Fig.  5.2-47)..  whereafter  tne 
expansion  part  of  the  wave  train  decreases  this  value.  The  shock  is  now 
displaced  faster  upstream,  although  the  post  shock  pressure  is  now 
decreasing  according  to  the  Imposed  wave  train. 

The  shock  reaches  its  minimum  location  at  t- 1  Orr/6  (Fig.  5.2-412), 
whereafter  it  starts  to  move  downstream.  At  t~(  1 5- 1  6)tt/6  it  has  again 
reached  its  steady  state  position  (Figs.  5.2-4.17  and  5.2-4.18).  which 
approximative^,  for  this  sample  case,  corresponds  to  the  maximum  shock 
position  for  the  first  cycle. 

The  procedure  is  then  repeated,  although  with  some  slight  differences  in 
the  time,  shock  location  and  pressure  level  as  the  flow  adapts  itself  to 
the  periodic  pressure  fluctuation  at  the  outlet  (these  plots  are  not  shown 
in  the  sequence  of  figures).  Then,  after  about  two  full  cycles,  the  flow  is 
periodic.  Such  a  periodic  cycle  is  shown  in  the  sequence  of  Figs.  5.2-4.20 
to  5.2-4.32.  The  minimum  shock  position  (Fiq.  5.2-4.23.24)  is  reached  at 
t*33.5n/6  (which  corresponds  to  a  time  lag  towards  the  exit  pressure  of 
M*0.8-tcycie)  and  the  maximum  (Figs.  5.2-4.29,30)  at  t*39.5n/6  (time  'ag 
towards  exit  pressure  At*0.3-tcyc]e). 

As  this  sample  case  falls  into  the  cathegory  tcycie/tsh~0(l)  according  to 
the  classification  in  section  5.2.1,  we  would  expect  the  minimum  shock 
position  to  appear  at 

t  =[time  for  disturbances  to  travel  from  the  exit  to  the  shock] 

+[tt/2  phase  lag  of  shock  position  towards  post  shock  pressure] 
♦[time  for  exit  pressure  to  reach  its  maximum  level] 
*TT/2+TT/2+TT/2=0.75-2TT=0.75tcycie 

and  the  maximum  at  1/2  period  later,  i.e.  at  t*0.25  tcycie.  it  is  seen  mat 
both  of  these  values  agree  well  with  the  numerically  determined  results. 
Again  it  should  be  pointed  out  that  the  fully  developed  periodic  solution 
has  a  maximum  shock  postion  xSimax«0.92  and  a  minimum  xs  mm*0.83, 
compared  to  the  steady  state  solution  of  xS|Steady=  1  °0. 
it  is  also  interesting  to  note  that  although  the  maximum  pressure  is  sucn 
that  it  would  nave  given  subsonic  flow  throughout  tne  nozzle  in  steady 
state,  the  shock  still  stays  well  behind  the  throat  (vthroot=0-5Cn  dirmg 


the  whole  period.  This  Is  clearly  a  fact  because  the  fluctuation  frequency 
is  in  the  same  order  of  magnitude  as  the  time  tSh 
Furthermore,  as  showed  in  section  5.2.1,  although  the  shock  is  in  the  same 
position  twice  during  its  oscillation,  the  pressure  distributions  at  these 
times  are  not  the  same.  This  Is  again  illustrated  by  the  average  shock 
position.  At  time  t=3 » tt/6  (Fig.5.2-4.2l)  the  shock  is  moving  forward  and 
at  t=37 n/6  (1/2  period  later,  Fig.  5.2-4.27)  it  is  moving  downstream.  The 
shock  location  is  approximatively  the  same,  but  the  unsteady  lift  and 
moment  introduced  downstream  of  the  shock  are  very  different. 

Most  of  the  difference  in  the  pressure  jump  across  the  shock  can,  in  this 
sample  case,  be  attributed  to  the  shape  of  the  imposed  wave  train,  and 
some  to  the  change  in  relative  shock  velocity  (see  section  5.2.1). 

Sample  Case  4. 

if  the  pressure  amplitude  is  increased  from  Pe,amp=005  to  De>amp=0.iO. 
with  all  the  other  parameters  the  same,  the  resulting  flow  through  the 
nozzle  will  be  completely  different.  We  will  here  only  discuss  the 
periodic  result,  beginning  at  the  time  t=5n  after  the  imposed  pressure 
fluctuation  at  the  exit  started  (Ftg.5.2-5. 1 ).  At  this  time,  the  exit 
pressure  is  decreasing  and  has  its  average  level.  The  shock  is  situated^ 
xs~0.83  at  an  intermediate  location  (Fig.  5.2-5. 3)  for  the  geometry  under 
consideration.  As  in  the  previous  sample  case,  the  steady  state  shock 
position  is  xSiSteady=!  .00. 

As  the  exit  pressure  decreases  further,  the  shock  moves  upstream  (Figs. 
5.2-5.4,5)  and  the  pressure  Jump  across  the  shock  decreases.  The  exit 
pressure  reaches  its  minimum  value  at  M1tt/2,  but  the  post  shock 
pressure  continues  to  decrease  because  of  the  time  lag  between  the  exit 
and  post  shock  pressures,  counteracted  by  the  increase  in  pressure  jump 
as  the  shock  moves  upstream,  i.e.  the  relative  velocity  the  snock  sees  is 
larger  than  in  a  steady  state  case  (Fig.  5. 2-5. 6).  The  pressure  just 
downstream  of  the  shock  does  therefore  not  decrease  as  much  as  would  oe 
expected  from  the  expansion  only.  This  can  be  the  explanation  for  the 
local  peak  in  post  shock  pressure.  It  is  also  probable  that  the  plateau  at 
point  "3"  in  the  figures  is  the  plateau  of  maximum  exit  pressure  that  has 
now  been  transported  to  be  located  just  downstream  of  the  shock. 


Simultaneously,  the  shock  has  moved  close  to  the  throat  ana  the  flow  is 
oarely  supersonic  upstream  of  it.  However,  tne  pressure  jump  is  still 
considerate.  Thereafter  tne  pressure  at  the  exit  increases  again,  ano  a 
compression -wave  moves  upstream  from  the  exit  (Figs.  5. 2-5. 7, 6.9' 
Simultaneously,  the  shock  has  now  moved  upstream  of  the  throat  and  it  is 
slowly  weakened  Py  the  decrease  in  flow  velocity  upstream  of  it  (as  the 
area  chances  the  relative  shock  velocity  will  become  subsonic),  as  well  as 
by  the  fact  that  the  post  shock  pressure  is  still  decreasing  (due  to  the 
finite  time  it  takes  for  the  disturbances  from  the  exit  to  propagate 
upstream)  it  is  aiso  smoothed  out  towards  a  compression  wave  isee  later 
discussion).  Between  the  shock  and  the  last  compression  wave  at  the  exit 
the  pressure  distribution  indicates  a  complicated  mixture  of  compression 
and  expansion  waves.  The  flow  immediately  downstream  of  the  shock 
builds  now  up  the  velocity  on  its  way  to  expand  towards  the  steady  state 
solution,  which  it  would  eventually  reach  if  the  exit  pressure  was  now  to 
be  kept  constant  at  the  average  value  (=level  of  Fig.  5.2-5. 9). 

Downstream  of  this  a  small  compression  and  another  expansion  lies  before 
the  large  compression  wave  at  the  outlet.  A  part  of  this  consists  of  tne 
upstream  moving  signals  being  reflected  at  the  shock,  a  part  by  the  phase 
lag  between  tne  response  of  the  shock  to  the  post  shock  pressure  changes, 
and  a  part  can  be  attributed  to  the  increase  in  pressure  jump  over  the 
upstream  moving  unsteady  shock. 

As  the  time  increases  (Figs.  5.2-5.10,11,12),  the  upstream  moving 
unsteady  shock  wave  (by  now  so  smoothed  out  that  It  can  be  considered  as 
a  compression  wave  instead)  approaches  the  Inlet  and  the  fiow  behind  it 
builds  up  towards  cnoking.  Downstream  of  this  build  up,  a  small 
compression  wave  is  forced  to  steepen,  probably  because  of  the  mam 
expansion  build  up  and  the  compression  build  up  at  the  exit.  The  final 
compression  steepens  fairly  soon  into  a  shock  wherefore,  at  t=39n/6, 
three  small  shocks  can  be  found  at  three  different  locations  in  the  nozzle. 
A  further  increase  in  time  pushes  the  first  shock  out  through  the  inlet 
(Figs.  5.2-5.13,14,15).  without  any  nottcable  reflection.  Simultaneously, 
the  main  expansion  is  built  up  and  the  two  last  snocks  merge  into  one 
large  instead,  in  tms  procedure  it  appears  that  tne  first  of  those  is 
moving  downstream  (pusned  by  the  main  expansion  build  up)  ano  tne 
second  is  moving  upstream  (pushed  by  the  imposed  pressure  at  the  outlet). 
Although  the  pressure  distripution  is  now  fairly  like  the  steady  state 
solution,  the  flow  is  still  not  sonic  at  the  throat.  The  blockage  appears 


somewhat  later,  at  t=44n/6,  Just  before  the  shock  is  again  pusheo 
upstream  of  the  throat  due  to  the  next  period  of  the  pressure  variation 
As  the  flow  is  now  periodic,  the  same  resulting  pressure  distributions 
will  repeat  themselves. 

In  the  present  sample  case  It  is  hard  to  define  an  “average  shock 
position",  but  it  is  still  of  interest  to  compare  two  times  during  the 
cycle.  As  the  reduced  frequency  is  the  same  as  in  sample  case  3,  the 
pressure  distributions  at  t=3ln/6  (Fig.  5.2-5.4)  and  t=37n/6  (Fig.  5^- 
5.10)  are  compared,  it  should  again  be  stated  that  these  two  positions  do 
not  correspond  to  tne  average  unsteady  exit  pressure. 

Tne  two  pressure  distributions  look  very  different,  which  again  (compare 
the  third  sample  case)  give  unsteady  lifts  and  moments  quite  different 
from  the  steady-state  values. 

it  should  also  be  pointed  out  that,  although  the  supersonic  shocks  are 
sharply  captured',  the  unsteady  shocks  In  the  subsonic  flow  are  spread  out 
a  bit  more. Two  reasons,  one  physical  and  one  numerical,  may  be  given. 
Firstly,  as  tne  snock  travels  upstream  the  relative  shock  velocity  may 
eventually  become  subsonic.  The  decrease  in  nach  number  then  will 
counteract  the  shock.  The  unsteady  upstream  moving  expansion  wave 
downstream  of  the  shock  (introduced  by  the  shock  moving  into  a  section 
with  increasing  area)  will  catch  up  with  the  shock,  and  it  will  start  to 
spread  out.  Secondly,  it  is  probable  that,  as  the  pre-shock  Mach  number  is 
subsonic,  the  numerical  dissipation  associated  with  a  supersonic  shock 
does  not  manifest  itself.  A  more  sophisticated  flux  limiter,  based  on  the 
gradients  instead  of  just  the  present  crude  way  of  determining  tne 
passage  over  the  sonic  line,  would  certainly  improve  the  subsonic  shock 
structure,  it  is  also  possible  that  the  van  Leer  flux  vector  splitting 
approach  may  improve  the  accuracy. 

Sample  Case  5. 

Obviously,  it  is  also  possible  to  drive  the  shock  from  its  steady  state 
position  into  the  subsonic  part  of  the  nozzle  while  keeping  the  exit 
pressure  fluctuation  amplitude  at  pe  amp=0.05,  as  in  the  third  sample  case. 

1  Indications  in  ref.  / 49,50/  are  that  the  unsteady  shocks  may  be  . 
captured  ever,  more  sharply  with  the  van  Leer  splitting  and  some  more 
sophisticated  flux  limiter. 


This  can  be  done  by  reducing  the  frequency  of  the  oscillation.  Here  we 
will  show  the  first  cycle  of  a  case  with  a  three  times  slower  fluctuation 
(Fig.  5.2-6.)  than  previously  (k=0.106  and  k=0.3l8  resp.).  Out  of  the 
sequence  four  time  levels  are  shown  (t=0,  t={  1  7/36]-2tt,  t=[26/36)  2n, 
t=[29/36]  2tt). 

As  for  sample  case  4  ti.«  change  in  exit  pressure  drives  the  shock 
upstream.  An  important  diffe.  ^nce  towards  the  previous  two  cases  :s 
though  that  the  post  shock  pressu.  5  distribution  is  now  much  flatter  (.or 
the  frequency  is  lower),  see  Figs.  5.2-C.2.3.  There  is  thus  no  second  shock 
oeing  built  up  before  the  first  has  dissapeared  through  the  inlet  due  to  tne 
compression  wave  from  the  exit.  However,  as  for  the  previous  sample 
case,  it  is  seen  that  the  flow  builds  up  to  become  supersonic  before  it  is 
choked  /Fig.  5.2-6. 3). 

The  flow  becomes  sonic  at  the  throat  (Fig.  5.2-6.4)  and  the  shock  appears 
as  a  combination  of  this  build  up  and  the  compression  from  the  exit  of  the 
nozzle.  The  fact  that  the  flow  becomes  supersonic  again  before  it  is 
choked  is  aiso  noted  by  looking  at  the  shock  position  (where  xs  is,  as 
before,  defined  as  the  value  where  the  Mach  number  within  the  numerical 
shock  structure  is  1)  in  Fig.  5.2-6. 5.  It  is  concluded  that  the  shock  is 
driven  forward  until  it  is  completely  subsonic  (at  t*35).  This  appears  at 
the  throat.  Thereafter  It  is  not  possible  to  compute  a  shock  location 
according  to  the  definition  that  Ms=1.00,  until  the  flow  once  again 
becomes  supersonic.  This  happens  at  t«4!.6  (compare  Fig.  5.2-6. 3),  before 
tne  flow  is  cnoxed  and  slightly  before  the  compression  wave  has  built  up 
to  a  shock.  Note  that  Ms=l.  first  appears  at  xs*0.72  (compare  Figs.  5.2- 
6.3,5),  i.e.  well  (Ax*0.22)  downstream  of  the  throat. 

As  the  exit  pressure  amplitude  is  larger  in  sample  case  4  than  in  3  and  5, 
the  shock  moves  faster  in  this  case.  For  the  two  cases  with  identical 
fluctuation  amplitude  the  average  shock  speed  is  approximative^  the 
same  (slightly  larger  for  the  highest  frequency).  This  can  be  expected 
from  the  fact  that  the  shock  velocity  in  a  shock  tube  increases  wirn 
increasing  pressure  ratio. 


Fig.  5.2-2a.  Computed  displacement  of  shock  position  in  time  for 
Adamson  /44/  second  sample  case,  k*0.035,  Pe.amp*00140 


RUN  NUMBER  -  96 


Fig.  5.2-2D.  Mach  number  distribution  for  steady  state  (at  t=4).  i  id,  ! 
x s=0  998)  as  well  as  approximate  minimum  (at  t=78.i2,  xs=0.8i0)  and 
maximum  (at  t=  129.52,  xs=l.215)  unsteady  shock  positions  during  the 
first  cycle.  Adamson’s  /44/  second  sample  case.  k=0.035,  pe  amp=00i40. 
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Fig.  5.2-3b.  Shock  velocity  for  k=0.035  (Adamson's  second 

sample  case,  all  parameters  as  In  Fig.  5.2-1).  Note:  Only  each  lOOth  time 
step  used  to  calculate  the  shock  velocity. 
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Fig.  5.2  -A. 27.  Unsteady  pressure  distribution  in  a 
Pe,amp=0.05.  All  other  parameters  as  in  Figs.  5.2-1. 
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Fig.  5.2-6. 1 .  Steady  state  pressure  distribution  in  a  nozzle  for 
Pe,amp=,")  05.  All  other  parameters  as  in  Figs.  5.2-1. 
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6.  Two  Dimensional  Unsteady  Results. 

Unsteady  two  dimensional  results  will  De  presentee  tor  a  channel  now 
wit*  a  1 0%  "Pump”  on  tne  upper  w3ii.  The  geometry  is  simi)ar  to  *r'p 
steady  state  case  presented  in  section  4.6.2,  apart  from  the  fact  that  'he 
"Pump"  has  new  a  sinusoidal,  instead  of  a  circular  arc,  form  This  change 
is  deliberately  introduced  in  order  not  to  create  any  numerical 
disturbances  at  the  corners  of  the  "bump"  in  cases  when  the  unsteady 
snock  propagates  thereover. 

Two  sample  cases  will  oe  presented.  Both  cases  consider  tne  same  steaov 
state  flow  field  as  in  section  4.6.2,  i.e.  pe=0.7369  (14-^=0.675).  The  first 
treats  an  exit  pressure  amplitude  of  Pe,amD=0- 12,  with  a  reduced  frequency 
(calculated  with  qr  as  in  the  previous  section)  of  k=0.396  (corresponding 
to  f=369Hz  for  a  profile  length  of  lr=0. Im).  The  second  case  has  a  iower 
amplitude,  pe  arrip=0.06,  with  the  same  reduced  frequency 
in  Doth  sample  cases  the  imposed  exit  pressure  is  the  same  over  the 
width  of  the  channel,  and  the  variation  of  pressure  starts,  as  in  the  quas’ 
one  dimensional  cases,  with  a  sinusoidal  increase. 


in  tne  first  sample  case  we  will  present,  starting  from  tne  beginning  of 
the  variation  of  exit  pressure,  about  1.5  periods.  The  steady  state  solution 
is  given  in  fig.  6-1,  and  the  unsteady  sequence  is  studied  in  Figs.  6-2  to 
6-31.  It  should  be  pointed  out  that  the  grid  used  for  this  computation 
consists  of  (14x60)  cells,  and  the  shock  is  captured  within  two  cells,  i.e. 
with  one  intermediate  computational  point  (see  for  example  Fig.  6-1). 

As  the  exit  pressure  increases,  the  lines  of  constant  Mach  number  close  to 
tne  outlet  are  pushed  upstream,  especially  at  tne  lower  waii  (Figs.  6-2  to 
6-4).  During  this  time  there  is  no  change  in  the  flow  up  to  x~0.~.  At 
t=ir/2  the  exit  pressure  has  reached  its  maximum  value  'Fig.  6-5', 
whereafter  it  decreases  (Fig.  6-6  to  6-8).  The  lines  of  constant  Mach 
number  close  to  the  exit  now  become  turned  in  the  other  direction 
(towards  the  upper  wall),  and  the  shock  starts  simultaneously  to  move 
upstream.  The  maximum  Mach  number  upstream  of  the  shock  decreas-rs 
raoidiv  as  the  shock  moves  upstream.  At  the  same  time  the  veioutv 
decreases  also  on  me  iower  waii  of  tne  nozzle,  ana  tne  isomacn  imes 


which  where  cent  towards  the  outlet  in  the  steady  state  (see  for  examoie 
m=0  70  in  Fig  6-1 )  are  no  dent  forward  instead  (M=0.70  in  Fig.  6-6). 

Ar  time  t=r>  (Fig.  6-9)  the  exit  pressure  has  agam  reached  its  steady 
state  level  out,  due  to  the  time  lag  between  the  pressure  at  the  exit  ano 
the  pressure  in  the  nozzle,  the  velocity  throughout  the  nozzle  continues  to 
decrease.  At  this  time  the  shock  is  pushed  into  the  subsonic  part  of  the 
nozzle  and  the  flow  is  completely  subsonic.  At  an  approximate  time  of 
t~9-2n/i6  the  velocity  has  reached  its  minimum  value  throughout  the 
nozzle,  and  the  flow  starts  finally  to  react  to  the  previous  decrease  in 
out  let  pressure  ne  tne  decrease  that  started  at  t=n/2,  Fig.  6-5).  Thus 
the  flow  velocity  begins  to  increase.  Before  this,  as  in  one  dimensional 
flow,  the  unsteady  shock  moves  upstream  in  the  subsonic  flow  domain 
(cigs.  6-9  to  6.10).  It  is  however  not,  in  contrast  to  the  quasi-one 
dimensional  results  from  section  5.2.2,  pushed  all  the  way  to  the  inlet, 
but  reaches  instead  a  minimum  position  and  is  then  moved  downstream, 
starting  at  t*9-2n/16. 

The  increase  in  flow  velocity  appears  over  the  whole  channel  width,  out 
obviously  witn  a  larger  value  at  the  bump.  As  the  imposed  decrease  in  rne 
pressure  wave  tram  approaches  the  bump,  the  flow  velocity  at  the  unpe- 
wall  will  again  reach  a  sonic  value  (t=l  1-2tt/!6,  Fig.  6-12),  and  the  snor 
steepens  again  (Fig.  6-13).  At  time  t=3n/2  (Fig.  6-13)  the  exit  pressure 
has  reached  its  minimum  level,  and  at  time  t~  1 3  2n/ 1 6  (Fig.  6-14)  tne 
shock  has  reached  its  steady  state  position,  however  with  a  lower  -re- 
shock  Mach  number  and  a  lower  pressure  jump  than  in  the  steed*  re¬ 
flow  field.  The  pressure  at  the  outlet  is  now  increasing,  but  the  ;  ;* 

still  moving  downstream  (Decause  of  the  finite  time  it  ta*-s  •  • 
disturbances  to  travel  upstream)  The  shock  at  the  upper  wa  1  •  - 

maximum  position  at  t*!9-2TT/l6,  with  a  maximum  pre  s^oo  ’  • 
significantly  higher  than  the  corresponding  steady  c*.r- 
maximum  position  is  reached  well  after  the  exit  re."  - 
reached  its  average  (=steady  state)  level  a: 

Simultaneously,  the  pressure  increase  at  the  cm’.  . 
and  this  compression  wave  is  slowlv  buiidma 
lower  wail  also  ana.  indeed,  over  rrv 
(C=l9-2n/i6,  Fig.  6-20).  Although  rrv?  --  . 

wall,  it  can  still  be  considered  to  •  ■ 

r  9  ■  lc  ( i- — ooo-r-t  /  \  <2  6_°'l  c. 


'earn  an 


s  ^1-  i  '  C- 


lower  wall  is  still  moving  uostream,  but  is  in  the  subsonic  field  of  the 
fiow.  At  t=2i -2n/i6  tne  snocK  at  tne  upper  wail  has  again  reached  us 
steady  state  location  (half  a  period  after  it  reached  the  same  position 
moving  downstream>.  The  pre-shock  Mach  number  is  now  much  higher  than 
the  corresponding  value  for  the  downstream  moving  shock  (Fig.  6-14),  and 
also  considerably  higher  than  the  steady  state  value  (Fig.  6-1).  The  rest  of 
the  flow  field  in  the  nozzle  is  also  very  different  between  these  three 
time  levels. 

The  shock  at  the  iower  wail  is  pushed  further  into  the  subsonic  part  of 
the  nozzie  ana  becomes  weaker  because  of  the  decrease  in  upstream  fiow 
velocity,  at  the  same  time  as  the  shock  at  the  upper  wall  also  is  moved 
upstream  (Figs.  6-23  to  6-25).  The  upper  part  of  the  shock  is  thereafter 
finally  again  pushed  into  the  subsonic  part  of  the  nozzle  (Fig.  6-26).  The 
shock  that  extended  over  the  whole  nozzle  width  at  time  t=20-2TT/l6  (Fig. 
6-21)  is  now  so  weak  that  it  is  only  noticed  as  a  shock  (unsteady, 
although  in  the  subsonic  flow)  close  to  the  upper  wall.  As  we  are  now 
considering  the  second  period  of  the  oscillation,  the  fiow  field  in  the 
nozzle  is  far  from  identical  to  the  flow  during  the  first  cycie  (compare 
Fig.  6-10,  6-26).  However,  as  the  shock  now  again  begins  to  move 
upstream  at  t=26-2rr/i6,  the  flow  field  is  slowly  built  up  towards  a 
periodic  solution.  Thus,  the  flow  during  the  second  cycle  is  fairly  close  to 
the  flow  field  of  the  first  cycle  (compare  Figs.  6-12,  6-28;  6-13,  6-29; 
and  6-14,  6-30),  and  the  same  pattern  of  shock  build  up  as  seen  during  the 
first  cycie  is  repeated.  In  fact,  the  flow  becomes  periodic  so  fast  that 
already  the  second  period  can  almost  be  considered  as  the  periodic 
solution.  This  can  oe  seen  oy  comparing  the  flow  field  at  t=25  2tt/  1 6  and 
one  period  later,  at  t=4i  -2tt/  i 6  (Figs.  6-26  and  6-31),  as  weii  as  tne 
corresponding  solution  at  later  periods  (not  shown  here). 


The  second  sample  case  treats  the  flow  in  the  same  nozzle,  with  the  same 
steady  state  solution  (iL^O.675.  pe=0.7369)  and  the  same  reduced 
frequency  of  the  imposed  exit  pressure  fluctuation.  The  difference 
towards  the  first  sample  case  is  the  amplitude  of  the  pressure  oscillation 
which  is  now  Peiamo=0.06,  towards  the  previous  Pe.amp^1?- 


««  **  V  V V 


The  steady  state  flow  field  is  given  in  Fig.  6-32,  and  a  few  time  steps  of 
the  first  unsteaoy  period  are  presented  in  Figs.  6-33  to  6-42. 
as  the  pressure  increases  at  the  exit  the  Isomach  lines  are  again  turned, 
as  in  the  first  sample  case,  upstream,  although  to  a  lesser  extent  as  the 
pressure  increase  is  smaller  (compare  Figs.  6-5  and  6-33).  As  the 
physical  time  for  the  disturbance  wave  train  to  travel  from  the  exit  tc 
the  shock  is  the  same  in  both  sample  cases,  the  shock  will  start  to  react 
at  the  same  time.  It  will  however  move  much  slower  in  the  case  with  the 
lower  pressure  amplitude.  This  is  confirmed  by  comparing,  for  example 
Figs.  6-6  ana  6-34. 

At  time  t=tr  the  exit  pressure  is  again  at  its  steady  state  value  (Fig.  6- 
35),  and  at  t*l0-2n/16  the  shock  at  the  upper  wall  has  reached  its 
minimum  position  (Fig.  6-36).  This  is  at  approximative^  the  same  time  as 
the  shock  had  reached  its  minimum  location  in  the  first  sample  case  (Fig. 
6-10).  There  are  a  fundamental  difference  though  between  the  high  and 
low  pressure  amplitude  cases.  As  was  concluded  for  pe  amp=0.l2  the  pre- 
snock  Mach  numoer  was  well  below  1.00  before  the  shock  started  to  move 
downstream  again  (Fig.  6.10),  whereas  for  Pe,amp=006  the  pre-shock  Mach 
number  is  still  supersonic,  and  the  supersonic  2one  around  the  location  of 
maximum  bumpthickness  is  clearly  manifested  (Fig.  6-36).  The  flow 
velocity  in  the  nozzle  starts  thereafter  to  increase  because  of  the  the 
decrease  In  exit  pressure  that  started  at  t=n/2.  The  shock  at  the  upper 
surface  of  the  bump  picks  up  in  strength  (Fig.  6-37  to  6-41),  and  reaches 
its  maximum  position  at  t*202t7/ 16,  sligthiy  after  the  time  tne  higher 
pressure  amplitude  shock  reached  its  maximum  location  (Fig.  6-20).  as 
expected,  for  the  lower  amplitude  case  the  compression  from  the  outlet 
can  not  build  itself  up  to  a  shock  wave  at  the  lower  wall,  wherefore  the 
flow  is  completely  subsonic  there  (Fig.  6-41).  This  is  in  contrast  to  the 
nigh  amplitude  case  where  the  flow  was  choked  during  a  small  part  of  the 
period  (Fig.  6-21). 

The  highest  pre-shock  Mach  number  (Fig.  6-41)  is  also  in  this  sample  case 
higner  than  the  steady  state  value,  although  not  as  high  as  in  tne  mgn 
amplitude  case  (Fig.  6-21),  and  the  supersonic  zone  reaches  out  to  *66% 
from  the  upper  wall  (see  M=1.00  in  Fig.  6.41),  compared  to  only  *44%  >n 
the  steady  state  solution. 


number  in  physical  plane 


Fig  6-1.  Sinusoidal  ”bump"-flow.  M-co=0.675,  p0=O.7369. 

Steady  state  solution  for  first  two  dimensional  sample  case. 


number  in  physical  plane 


Mg.  6-2.  Sinusoidal  ”bump"-flow.  M-oo'0.675,  pe=0.7369. 

Unsteady  solution  for  first  two  dimensional  sample  case. 


number  in  physical  plane 


Fig.  6-3.  Sinusoidal  "bump“-flow.  M.^0.675,  pe=0.7369. 

Unsteady  solution  for  first  two  dimensional  sample  case. 


number  in  physical  plane 


Fig.  6-4.  Sinusoidal  "bump"-flow.  M-o^O.675,  pe=0.7369. 

Unsteady  solution  for  first  two  dimensional  sample  case. 
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Fig.  6-6.  Sinusoidal  ”bump“-flow.  M-oo'0.675,  pe=0.7369. 

Unsteady  solution  for  first  two  dimensional  sample  case. 
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Fig.  6-7.  Sinusoidal  "bump"-flow.  M-^0.675,  pe=0.7369. 

Unsteady  solution  for  first  two  dimensional  sample  case. 


Fig.  6-8.  Sinusoidal  "bump"- flow.  11.^=0.675.  oe=0.7369. 

Unsteady  solution  for  first  two  dimensional  sample  case. 
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Fig.  6-9.  Sinusoidal  "bump'’-flow.  11-^=0.675,  pe=0.7369. 

Unsteady  solution  for  first  two  dimensional  sample 
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Fig.  6-10.  Sinusoidal  "bumpB-f!ow.  M.^0.675,  pe=0.7369. 

Unsteady  solution  for  first  two  dimensional  sample  case. 


Fig.  6-11.  Sinusoidal  "bump"-flow.  M-^0.675,  pe=0.7369. 

Unsteady  solution  for  first  two  dimensional  sample 


ach  number  in  physical  plane 


Fig.  6-12.  Sinusoidal  "bump”-flow.  M.^0.675,  pe=0.7369. 

Unsteady  solution  for  first  two  dimensional  sample  case. 


ach  number  in  physical  plane 


Fig.  6-13.  Sinusoidal  "bumpWlow.  M.^0.675,  pe=0.7369. 

Unsteady  solution  for  first  two  dimensional  sample  case. 
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Fig.  6-14.  Sinusoidal  Mbump”-flow.  M.^0.675,  pe-0.7369. 

Unsteady  solution  for  first  two  dimensional  sample  case. 
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Fig.  6-15.  Sinusoidal  "bump"-flow.  M-^0.675,  pe=0.7369. 

Unsteady  solution  for  first  two  dimensional  sample  case. 
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Fig.  6-16.  Sinusoidal  "bumpM-flow.  M-^0.675,  pe=0.7369. 

Unsteady  solution  for  first  two  dimensional  sample  case. 
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Fig.  6-17.  Sinusoidal  “bumpM-flow.  M.^0.675,  pe=0.7369. 

Unsteady  solution  for  first  two  dimensional  sample  case. 
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Fig.  6-19.  Sinusoidal  "bump"-flow.  M.^0.675,  pe=0.7369. 

Unsteady  solution  for  first  two  dimensional  sample 
Pft.amp=Q  i2,  k=0.369. _ ; _ i__ _ 
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Fig.  6-20.  Sinusoidal  ”bump"-flow.  M.^0.675,  pe=0.7369. 

Unsteady  solution  for  first  two  dimensional  sample  case. 
Pe,amp=0- 1 2,  k=0.369. 


number  in  physical  plane 


'  Fig.  6-22.  Sinusoidal  "bump"-flow.  M_«,=0.675(  pe=0.7369. 

unsteady  solution  for  first  two  dimensional  sample  case.  . 

Pe,amp=Q' ' ->  k=0^69. 
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Fig.  6-24.  Sinusoidal  "bump'-flow.  M-^0.675,  pe=0.7369. 

Unsteady  solution  for  first  two  dimensional  sample  case, 
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Fig.  6-?6.  Sinusoidal  "bump"-flow.  M.^0.675,  pe=0.7369. 

Unsteady  solution  for  first  two  dimensional  sample  case 


Fig.  6-27.  Sinusoidal  "bump" -flow.  M.^0.675,  pe=0.7369. 

Unsteady  solution  for  first  two  dimensional  sample  case. 


Fig.  6-28.  Sinusoidal  "Pump"-flow.  M.^0.675,  pe-0.?369. 

Unsteady  solution  for  first  two  dimensional  sample  case. 


Fig.  6-29.  Sinusoidal  "bump"-flow.  M.^0.675,  pe=0.7369. 

Unsteady  solution  for  first  two  dimensional  sample  case. 
Da  airin— 0. 1  2,  k=0.369. 


Sinusoidal  "bump'-flow.  fi.^0.675,  pe=0.7? 
Unsteady  solution  for  first  two  dimensional 
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Fig.  6-31.  Sinusoidal  Mbump"-flow.  M-^0.675,  pe=0.7369. 

Unsteady  solution  for  first  two  dimensional  sample  case. 

P8.amp=0-'2,  k=0.369. 
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Fig.  6-32.  Sinusoidal  "bump"-f)ow.  M-^0.675,  pe=0.7369. 

Steady  state  solution  for  second  two  dimensional  sample 
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Sinusoidal  "bump"-flow.  11-^=0.675,  pe=0.7369. 

Unsteady  solution  for  second  two  dimensional  sample  case 


Mach  number  in  physical  plane 


Sinusoidal  “bump"-flow.  M.^0.675,  pe=0.7369. 

Unsteady  solution  for  second  two  dimensional  sample  case 
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Fig.  6-39.  Sinusoidal  "bump’-flow.  11-^=0.675,  pe=0.7369. 

Unsteady  solution  for  second  two  dimensional  sample  case. 
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Fig.  6-41.  Sinusoidal  "bump'-flow.  11-^=0.675,  pe=0.7369. 

Unsteady  solution  for  second  two  dimensional  sample  case. 
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Fig.  6-42.  Sinusoidal  "bumpMlow.  M-^0.675,  pe=0.7369. 

Unsteady  solution  for  second  two  dimensional  sample  case. 


~  Transformation  From  the  Physical  Into  the  Computational 
Plane  of  Reference. 


Tr.e  general  transformation  from  the  (x,yft)-plane  to  the  (^^jj-plane  will 
Oe  presented  «n  this  section.  The  transformation  in  r  is  straightforward, 
as  r=t.  and  oniv  (x.y)  to  C$,i))  is  discussed. 

Tne  transformation  oresented  snouia  oe  user  rriendiv  in  as  far  as  it 

/ 

snouid  allow  for  simple  changes  of  cascade  geometry,  it  snouio  also  allow 
for  a  Place  vibration,  presently  in  circumferential  direction  (as  seen  m 
the  rotating  turbomachine,  y  direction  in  the  Cartesian  coordinate  system 
?r,  the  present  study),  and  easily  be  extended  to  accomodate  any  genera! 
blade  vibration.  In  order  to  treat  cases  of  propagation  of  disturbances 
through  a  cascade,  the  transformation  should  be  able  to  handle  several 
blade  passages  of  a  general  profile.  Furthermore,  the  mapping  should 
ailow  for  the  foliowing: 

The  inlet  and  outlet  planes  are  fixed  in  both  the  computational  ana 
physical  planes  during  a  generalized  blade  vibration. 

The  leading  edge  position  shall  be  constant  in  a  predefined  axial  and 
circumferential  direction  in  the  computational  plane. 

idem  for  the  trailing  edge. 

Computational  points  should  be  clustered  in  the  leading  and  trailing 
edge  regions,  and  close  to  the  blade  surfaces. 

Tne  computational  points  towards  tne  inlet  and  outlet  flow 
boundaries  in  the  physical  plane  should  be  stretched  away  from  the 
profile. 

The  mapping  should  be  fast  as  It  has  to  be  performed  for  every  time 
step  of  the  calculation,  both  for  the  calculation  with  vibrating  blades  and 
for  the  calculation  with  stationary  blades  as  in  both  cases  the  slipstream 
from  the  traillnq  edge  will  move  in  time. 


For  the  above  objectives,  it  was  decided  to  introduce  the  blade  geometry 
as  coordinate  points,  and  to  adopt  a  numerical  grid  of  "H"-form  (Fig  7.1). 
This  sort  of  grid  is  acceptable  for  capturing  normal  shock  waves  (wmcn 
is  the  scope  of  the  present  study).  Obviously,  a  C-gnd  or  an  O-grid  would 
be  better  from  the  point  of  view  of  accurately  resolving  the  flow  around 
the  leading  edge.  However,  the  creation  of  such  grids  for  unsteady  cascade 
flows  is  a  major  undertaking  in  itself.  It  is  instead  proposed  to  use  ar.  H- 
mesh  for  global  resolution,  and  in  the  future  overlap  this  global  H-mesh 


with  a  local,  around  and  fixed  with  the  airfoil,  C-mesh.  This  Kind  of 
approach  has  oeen  shown  by  Steger  /55/  to  be  feasible  for  steady  state 
riows  (Euler  solver)  and  for  small  unsteady  perturbation  potential  flow  py 
Verdon  /56/. 

The  transformation  includes  clustering  of  points,  for  better  resolution  in 
regions  of  high  flow  gradients,  in  several  locations,  presently  at: 

downstream  of  inlet  boundary  and/or  upstream  of  leading  edge  plane, 
downstream  of  leading  edge  and/or  upstream  of  trailing  edge  plane, 
downstream  of  trailing  edge  plane  and/or  upstream  of  outlet 
boundary. 

The  physical  domain  of  interest,  limited  by  the  inlet  and  outlet  boundaries 
and  the  upper  and  lower  blades  (Fig.  7.1),  is  then  separated,  in  ax'al 
direction  (=x),  into  three  different  regions,  each  extending  over  the 
number  of  blade  passages  choosen  in  the  circumferential  direction  (=y): 
Region  I:  From  inlet  to  leading  edge  plane  ($<0.). 

Region  2:  Blade  passages  (O.sgsl.). 

Region  3:  From  trailing  edge  plane  to  downstream  boundary  ($>i.). 

Each  of  these  regions  are  mapped  into  a  rectangular  computational  piane 
of  reference  with  the  following  transformations. 

The  mapping  in  axial  direction  (x)  is  given  by  eq.  7.1  in  regions  1  and 

2/1,2/: 

^(x,y,t)=^start+^end"^start^  (®x+^-®x^  ln(Aj/A2)/ln(A2/A^))  (/.i/ 

where: 

=px  "  2ax  +  (2ax+l)(x-xstart)/(xend-xstar{) 

A^  =Px  +  2<Xx  -  (2ax+  1  )(x-xstart)Axend"*start) 

=Px  +  1 
A^  =px  -  1 

a  =  variable  to  indicate  where  the  mesh  should  be  refined. 

=  0.0:  refinement  only  towards  the  end  of  the  region. 

=  0.5;  equal  refinement  both  towards  the  end  and  the  beginning, 
p  =  stretching  parameter  (p>l.) 

index  "x":  parameters  in  x-dlrection  (a/p  are  different  for  each 
region). 

index  "start" :  x  and  $  value  at  the  beginning  of  eacn  region, 
index  "end"  :  x  and  S  value  at  the  end  of  each  region. 


in  region  3  tnis  mapping  is  given  by  eq.  7.2: 


-$(x,y,t)=-*3tdrl+{$en(j-$startH- 1  +otx+(  1  -ax)ln(A9/AJ0)/1n(A3/A4)} 

(7.2) 

where: 


Ag  =  Px“2Ox+(20<x  +  I  )(X~XerNj)/(Xstap|-Xe|fHj) 

A10  =  Px*2o<x-(2ax+  1  )(X"X0n<j)/(x$tart~xen<p 


The  reason  for  this  change  in  the  mapping  between  regions  2  and  3  is  to 
have  the  possibility  to  cluster  points  in  the  beginning  of  zone  3,  instead 
of  towards  the  end  as  in  region  l. 

In  the  circumferential  direction  (y)  the  mapping  Is  given  by  eq.  7.3: 


n(x,y,t)=nst8rt+{nend-f|startKay+(  1  -«y)  ln(B, /B2>/ln(B3/B^)) 

(7.3) 


where: 


B1  =Py  -  2o<y  +  (2ay+l  )(y~ystart)/(yen<rystart) 

=Py  +  2(xy  -  (2ay+ l)(y-yS(art)/(yen(j-ystart^ 

-  B3  =py  ♦  1 

B«  =^y  '  1 

a  =idem  as  for  axial  direction. 

P  =idem  as  for  axial  direction, 
index  y :  parameters  in  y-direction  . 

index  "start”:  y  and  q  value  at  the  beginning  of  each  region  (i.e.  tne 
blade  surfaces). 

index  "end":  y  and  r\  value  at  the  end  of  each  region  (i.e.  the  blade 
surfaces). 


Here,  the  computational  grid  stretches  from  h=0.  (for  the  upper  surface  of 
blade  1)  till  r)=N-1  (for  the  lower  surface  of  blade  N).  Thus,  for  blade 
passage  "  1 ",  n$tart=0-  whereas  nend= 1  • 

This  transformation  maps  the  physical  plane  into  the  computational  plane 
(i.e.  (x,y)  to  (Ul»-  The  Inverse  of  the  mapping  ((^q)  to  (x,y))  is  easily 
computed  to  be: 


-  in  axial  direction,  regions  1  and  2: 

xi,j=*staft+(xend-><start)  f +  2ttxA7]/[A7(2<Xx+  I )]  (7.4) 

where 

i  :  computational  point  in  S,  direction 

j  :  computational  point  in  t|  direction 

A5  =  ^i,j-Sstart)A$end-£start) 

A6  =  (As-ax)/(l-ax) 

A?  =(px+l)A6  +  (px-l)A6 

a8  =CPx^)a6-Ox-i)A6 

-  in  axial  direction,  region  3: 

xi,j=xeod+^xstart_xend^PxAg  +  20(xA7l/[A7(2ocx+ 1 )]  (7.5) 

where  now  A5  is  different  from  regions  1  and  2: 

a5  =l~($i,j~$start)/^send"$start) 

-  in  circumferential  direction,  the  inverse  of  eq.  (7.3)  is: 

yi,j=ystart+^yend-ystart^  fPy  +  2o<yB7]/[B7(2o<y+ 1 )]  (7.6) 

where 

=  ^i.j^start^^end'Hstart^ 

Bg  =  (B5-ay)/(  1  -<Xy) 

B?  *  (py+  1  )Bg  +  (Py-  I  )Bg 

B,  =  (py+  I  )Bg  -  <Py"  1  )Bg 

An  example  of  the  mapping  here  presented  is  given  for  a  symmetrical 
Joukowski  profile  at  60°  stagger  angle  (from  axial)  in  Fig.  7.2.  Both  the 
entire  domain  of  consideration  (Fig.  7.2a)  as  well  as  blown  up  versions  in 
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different  part  of  the  flow  field  are  shown  for  the  physical  plane,  whereas 
a  general  view  of  me  computational  plane  is  presented  (Fig.  7.3) 
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F«g-  7.1  :  Physical  domain  of  reference  with  H-grid. 

Fig.  7.2  :6nd  of  jouxowski  profile  (15%  thickness)  at  60°  stagger  angle 
in  tne  physical  plane. 

Fig.  7.3:  Grid  of  Joukowski  profile  (15%  thickness)  at  60*  stagger  angle 
in  the  computational  plane. 

These  figures  can  be  seen  on  the  computer  terminal  by  running  the 
program  on  an  IBM- AT.  However,  no  hardcopy  device  exists  for  this  type  of 
application  in  the  "personal  Computer  Lab"  in  the  "Computer  Center”  at  the 
NPS.  Tne  figures  can  thus  not  be  shown  here. 


8  Numerical  Determination  of  Metric  Terms. 

The  metric  terms  are  calculated  as  centered  differences  (see  Fig.  7-i  > 

For  calculation  of  the  Jacobian  determinant  in  the  point  (j,i)  we  need 
values  from  the  surrounding  points. 

Dj,i=l /{WWJj.i 

with 

=  txj.i+rxj.i-i^ftj.i+r^j.i-i) 

Wj.i  =  tyj  .!♦  i  -yj  ,i- 1  >/Uj  .«♦  i  -^j  fi- 1 ) 

=  {Xj+j  ,pXj_j  ,j}/{T)j+j  rHj-1  ,i) 

fyiA  ,i =  Cyj+  i  ,ryj- 1  \  ,rnj- 1 ,0 

Furthermore,  for  the  “MUSCL'-approach,  the  derivatives  (xn,yn)  in  the 
points  (j, 1-1/2)  and  (j, M/2),  as  well  as  the  derivatives  (x^y*)  in  the 
points  (J+ 1/2.1)  and  (J- 1/2,1)  are  needed. 

These  are  evaluated  as: 

<X|j)j  .1+1/2  =  tXj+  j  j* 1  /2~Xj-  j  ,j*  j  /2)/(2Arj} 

=  (xj+ 1  ,i+  PXjh  ,rxj- 1  >  j-xj-i  (|}/{4Ar|) 

.1*1/2  =  Wj+lti+l/2-yj-J,1+1/2V{2All) 

=  fyj*  i  ,i+ 1 +yj+ 1  .ryj- 1  ,i+ 1  _yj- 1  .iJ/mati) 

(x^j  ,1-1/2  =  .i-1/2~xj-1  ,t-1/2>/t2ATl> 

-  txj^  I  f|*Xj+  j  j-  I  -Xj-  j  >rXj-  j  (j- 1  }/{4Ai]} 

fyr|}j,i-l/2  =  {yj+U-l/2‘yj-l.1-!/2V{2Ar|) 

=  (yj+ 1  ,ryj+ 1  ,i- 1  -yj- 1  ,ryj-  \  ,i- 1  >/(^An> 

m+1/2'.l  =  txj+  j  /2,|+ 1  -Xj+  j  /2,i- 1  >/  (2A$) 

=  (Xj+ 1 ,1+  i +xjii+ 1  *Xj+ 1  tj- 1  -Xjj- 1  )/{4A$) 

ty?}j*1/2,i  =  {yj+1/2.1+ryj+1/21l-lV(2A^) 

=  (yj+i>i+ryjti+ryj+i1i-ryj.i-iV(4^} 

^j-t/2,i  *  txj.  j  /2  ,t+ 1  ~xj- 1 /2 ,1- 1 )/  (2A^} 

=  (Xj  ,j+  j +Xj_  j  ,♦  j  -Xj  j.  j  -Xj.  |  |  )/{4AU 

{y*}j- 1/2.1  =  Cyj- 1  /2.1+- 1  -yj- 1  /2.I- 1  >/C2A^} 

-  (yj.i+ryj-u+ryj.i-ryj-u-iJ/t4^} 


l*U*rf«.ce 
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Fig.  8-1.  Grid  points  and  control  volume  for 
metric  terms. 


the  evaluation  of  the 


9.  Recommendations  for  Improvments. 

!n  the  following  a  few  suggestions  for  improving  the  accuracy  of  the 
results  obtained  by  the  present  method  are  given. 

!n  two  dimensions  the  crude  way  of  determining  the  departure  point 
of  the  characteristic  (u-a)  at  the  inlet  and  (u+a)  at  the  outlet  certainly 
introduces  some  numerical  reflection  (section  3.3.2).  In  ref.  /47/  an 
iterative  procedure  for  finding  this  departure  point  has  been  given.  This 
method  could  be  tried  out,  still  in  connection  with  the  "radiative' 
boundary  condition  at  the  inlet,  in  order  to  compare  the  numerical 
reflections  with  the  present  implementation. 

For  an  eventual  improvment  of  the  sharp  capture  of  unsteady  shocks 
it  would  be  of  interest  to  implement  the  flux  vector  splitting  in  van 
Leer's  approach,  and  to  compare  this  with  the  Steger/Warming  splitting, 
indications  in  the  literature  /50/  indicates  that  this  approach  may  give 
sharper  shocks  than  the  present  implementation. 

in  the  same  direction  it  would  be  of  interest  to  Iook  into  the 
possibility,  as  explained  by  Allmaras  and  Giles  in  /50/,  of  using  the  flow 
gradients  within  one  cell  only  for  obtaining  a  second  order  spatial 
accuracy,  instead  of,  as  present,  employing  two  grid  points  for  this 
(section  3.1 ). 

This  method  will  also  make  .it  possible  to  keep  a  second  order 
spatial  accuracy  at  the  blades,  without  the  need  for  extrapolating  any 
flow  variables.  Thus  the  accuracy  of  the  stagnation  enthalpy  at  the  blades 
wouid  improve. 

The  accuracy  of  the  stagnation  enthaply  will  also  improve  if  grid 
points  are  clustered  close  to  the  walls. 

For  a  better  accuracy  of  the  total  enthalpy  at  the  walls,  a  modified 
splitting  of  the  enthalpy  has  been  proposed  in  /59/.  The  modification 
should  be  small  and  fairly  easy  to  implement.  According  to  reference 
/59/,  the  improvement  of  accuracy  can  be  considerable. 

Quasi-one  dimensional  test  cases  by  Adamson  et  al.  (/39/-/4S/). 
with  shock  oscillations  so  large  that  the  shock  moves  upstream  of  the 
throat,  could  validate  such  numerical  solutions.  Comparisons  of  this  kind 
would  also  serve  to  better  understand  the  unsteady  flow  phenomena. 


In  order  to  validate  the  two  dimensional  unsteady  results,- 
comparisons  with  detailed  experiments  could  be  made.  Eventually  couio 
the  numerous  experimental  results  by  Sajben  et  a)  serve  as  test,  cases. 

An  interesting  direct  application  of  the  method  would  be  to 
consider  the  overlapping  part  of  a  transonic  compressor  cascade  as  3  two 
dimensional  curved  channel  and  to  investigate  how  large  an  unsteady  back 
pressure  fluctuation  must  be  in  order  to  significantly  move  the  blade 
passage  shock  away  from  its  steady  state  position. 
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